Classifications

• • G—PHYSICS
  • G06—COMPUTING; CALCULATING OR COUNTING
  • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
  • G06T11/00—2D [Two Dimensional] image generation
  • G06T11/003—Reconstruction from projections, e.g. tomography
  • G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
• • G—PHYSICS
  • G06—COMPUTING; CALCULATING OR COUNTING
  • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
  • G06T11/00—2D [Two Dimensional] image generation
  • G06T11/003—Reconstruction from projections, e.g. tomography
  • G06T11/008—Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
• • G—PHYSICS
  • G06—COMPUTING; CALCULATING OR COUNTING
  • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
  • G06T2211/00—Image generation
  • G06T2211/40—Computed tomography
  • G06T2211/421—Filtered back projection [FBP]
• • G—PHYSICS
  • G06—COMPUTING; CALCULATING OR COUNTING
  • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
  • G06T2211/00—Image generation
  • G06T2211/40—Computed tomography
  • G06T2211/424—Iterative

Definitions

• the invention relates to the field of image reconstruction.
• Major deficiencies associated with FBP include the need for a large number of projections to achieve limited quantitative accuracy.
• the number of projections is typically counted in the hundreds or thousands but the projections are not used as efficiently as they could be.
• IR IR
• the benefits of IR techniques include their ability to reduce reconstruction errors when following FBP.
• Finding the optimal object reconstruction requires the minimization of the objective function, which is the sum of all reconstruction errors. Finding the optimal reconstruction requires operations on each voxel. Near the minimum error the objective function can typically be represented by the Hessian matrix, describing the second derivative of the objective function.
• Multi-grid variations of these algorithms may help, but ultimately still fail because of the size of the Hessians involved with fine grids.
• Multi-grid resolution here refers to the use of progressively finer resolution as iterations are performed.
• characterizing for example, object deformation, beam hardening and scatter in the case of x-ray imaging, or field distortions in the case of magnetic resonance imaging (MRI), may also increase measurement data inconsistencies.
• MRI magnetic resonance imaging
• a system comprises: a non- transitory data storage for storing projection space data, the projection space data in a density domain for an object under observation, and including one or more input projection pixels and one or more predicted projection space pixels; and an image reconstructor computer having at least one processor, the at least one processor operable to: receive the projection space data in the density domain from the non-transitory data storage; compute one or more measured transformed pixels in a transformed domain using the one or more input projection pixels in the density domain and input transformation functions; compute one or more predicted transformed pixels in the transformed domain using the one or more predicted projection space pixels in the density domain and reference transformation functions; compute first pixel innovation result data in the transformed domain using a difference between the one or more measured transformed pixels and the one or more predicted transformed pixels; compute a pixel-by-pixel innovation scaling matrix using inverse slopes of the input transformation functions and inverse slopes of the reference transformation functions for corresponding input and reference pixel values; compute second pixel innovation result data using a pixel-by-pixel product of the
• a non-transitory data storage device stores software code executable by a computer having one or more processors, the software code to: receive projection space data in a density domain for an object under observation, the projection space data including one or more input projection pixels and one or more predicted projection space pixels; compute one or more measured transformed pixels in a transformed domain using the one or more input projection pixels in the density domain and input transformation functions; compute one or more predicted transformed pixels in the transformed domain using the one or more predicted projection space pixels in the density domain and reference transformation functions; compute first pixel innovation result data in the transformed domain using a difference between the one or more measured transformed pixels and the one or more predicted transformed pixels; compute a pixel-by-pixel innovation scaling matrix using inverse slopes of the input transformation functions and inverse slopes of the reference transformation functions for corresponding input and reference pixel values; compute second pixel innovation result data using a pixel-by-pixel product of the first pixel innovation result data and corresponding elements of the pixel-by-pixel innovation scaling matrix; compute preliminary transformed
• voxel-by-voxel update scaling matrix wherein at least one voxel of the transformed object voxel density update estimate is associated with an element of the voxel- by-voxel update scaling matrix; add the transformed object voxel density update estimate to a corresponding transformed preceding voxel data estimate to obtain a transformed density estimate; and reconstruct an object space image representing the object under observation using the transformed density estimate.
• a method for image reconstruction performed by an image reconstructor computer having at least one processor comprises:
• the projection space data including one or more input projection pixels and one or more predicted projection space pixels; computing one or more measured transformed pixels in a transformed domain using the one or more input projection pixels in the density domain and input transformation functions; computing one or more predicted transformed pixels in the transformed domain using the one or more predicted projection space pixels in the density domain and reference transformation functions; computing first pixel innovation result data in the transformed domain using a difference between the one or more measured transformed pixels and the one or more predicted transformed pixels; computing a pixel-by-pixel innovation scaling matrix using inverse slopes of the input transformation functions and inverse slopes of the reference transformation functions for corresponding input and reference pixel values;
• computing second pixel innovation result data using a pixel-by-pixel product of the first pixel innovation result data and corresponding elements of the pixel-by-pixel innovation scaling matrix; computing preliminary transformed object update data using a tomographic
• FIG. 1 depicts an imaging system according to an exemplary embodiment of the invention
• FIGS. 2A and 2B depict example data input, reference prediction, and inner feedback loop 165 of enhanced high efficiency computer tomography (CT) with optimized recursions with reconstructed object density according to an exemplary embodiment at iteration i;
• CT computer tomography
• FIG. 3 depicts an example x-ray projection during contrast dye injection
• FIG. 4 depicts a reconstruction of a beating hydraulic coronary tree model from six projections according to an exemplary embodiment of the Levenb erg-Mar quardt process applied to 3D cone beam reconstruction;
• FIG. 5 depicts an example influence function
• FIGS. 6A and 6B depict example reconstructions of a head phantom from 200 projections (FIG. 6A) and five projections (FIG. 6B).
• FIG. 1 depicts an imaging system according to an exemplary embodiment of the invention.
• the imaging system may include: a controller 135; an interface 140 (e.g. a graphic user interface, keyboard, joystick, etc.) for signal communication with an investigator 170 (e.g., a computer) and synchronization with the controller 135; transmitted s) 115 for generating and emitting excitation energy 120 in response to a control signal from the controller 135; and detector(s) 145 configured to generate measured pixel data 155.
• the measured pixel data 155 may be in a digital format stored in a data storage device 180.
• Measured pixel data 155 containing information encoding the internal structure of an object 125 may be transformed using image reconstructor 165 (e.g., a computer) into a reconstructed image data 261 to be visualized on an output device 175, for example.
• the measured pixel data 155 may be from an experimental acquisition 110, a simulation 185 (e.g., on a computer), and/or a data storage device 180 containing recorded projected pixel data from an earlier experiment or simulation, for example.
• the experimental acquisition 110, data storage device 180, and simulation 185 may remotely provide the measured pixel data 155 to the image reconstructor 165 directly or via a network, for example, a local area network (LAN), a wide area network (WAN), the Internet, etc.
• LAN local area network
• WAN wide area network
• the Internet etc.
• Measured pixel data 155 may include the result of delivering any form of excitation energy into an object under observation 125 and thus may include data from, for example, electron microscopy, magnetic resonance (MR), positron emission tomography (PET), single positron emission computed tomography (SPECT), ultrasound, fluorescence, multi-photon microscopy (MPM), optical coherence tomography
• MR magnetic resonance
• PET positron emission tomography
• SPECT single positron emission computed tomography
• MPM multi-photon microscopy
• Data storage device 180 may be, for example, non-transitory memory, CD-ROM, DVD, magneto-optical (MO) disk, hard disk, floppy disk, zip-disk, flash-drive, cloud data storage, etc.
• Measured pixel data 155 from experimental acquisition 110 and/or simulation 185 may be stored on data storage device 180.
• the imaging system may also include: a translating or rotating table 130 configured for receiving an object 125 thereon and operable to translate or rotate in relation to the transmitted s) 115 and the detector(s) 145; and an image reconstructor 165 coupled electrically, or via an optional data transporter 150, such as the Internet, other networks, or a data bus, to the detector(s) 145 and the controller 135.
• a data bus may be, for example, a subsystem of electrical wires that transfers data between computer components inside a computer or between computers. Although the data flow between connected components is uni-directional, the communication between connected components may be bi-directional.
• Excitation energy 120 may be a form of radiation energy such as, for example, X-ray energy, electromagnetic (EM) energy, optical energy, infra-red (IR) energy, particle energy (e.g., electron, neutron, atom beams), vibration energy, such as ultrasound, etc.
• Excitation energy 120 may be irradiated onto an object 125, which may be, for example, a phantom, a human patient, a specimen, or any combination thereof.
• the excitation energy 120 may be emitted by transmitter(s) 115 of the corresponding energy type.
• the excitation energy 120 may propagate through the object 125 and a portion may be received by the appropriate detector(s) 145.
• the detector(s) 145 may convert received energy into measurable electrical signals that may be further convert the measured electrical signals into projected pixel data in a digital format.
• Controller 135 may be a circuit connecting the transmitter(s) 115 and detector(s) 145 and may send control signals to the transmitted s) 115 and detector(s) 145 to synchronize the transmission of the excitation energy 120 and the operation of the detector(s) 145.
• the circuit may be analog, digital, or mixed-signal.
• Controller 135 may also be a computer having one more processors and one or more memories to send control signals to the transmitted s) 115 and detector(s) 145.
• the image reconstructor 165 may be responsive to the controller 135 and receptive of the measured pixel data 155 to reconstruct an image of the object 125 via a method according to some embodiments of the invention to produce a high fidelity image of the object with high computation efficiency.
• the image reconstructor 165 may be, for example, a computer having one or more data storage devices storing software code to operate the computer according to the exemplary embodiments of the invention.
• the computer may include one or more processors and non-transitory computer-readable medium to read and execute the software code stored on the one or more data storage devices.
• the image reconstructor 165 may be, for example, a computer having one or more data storage devices storing software code to operate the computer according to the exemplary embodiments of the invention.
• the computer may include one or more processors and non-transitory computer-readable medium to read and execute the software code stored on the one or more data storage devices.
• the reconstructor 165 may include one or more program storage and execution devices operable according to the exemplary embodiments of the invention.
• the image reconstructor 165 may produce the reconstructed image data 261.
• the image reconstructor 165 may receive the measured pixel data 155 and may process the measured pixel data 155 by, for example, nonlinear transformation 200 to generate cumulative object density projection data 201, further computing measured transformed pixels 206, the nonlinear transformations of projection data 201, and computing their excess over predicted transformed pixel data 211, to generate innovations, to produce reconstructed image data 261, and may further generate, for example, predicted projection data 276 in Figure 2.
• An output device 175 may receive one or more of the reconstructed image data 261, reconstructed image error data 160, or reconstruction image error data 216.
• Output device 175 may be a visualization device or a data storage device, for example.
• a visualization device may be, for example, a display device or a printing device.
• Example display devices may include, for example, a cathode ray tube (CRT), a light-emitting diode (LED) display, a liquid crystal display (LCD), a digital light projection (DLP) monitor, a vacuum fl orescent display (VFDs), a surface- conduction electron-emitter display (SED), a field emission display (FEDs), a liquid crystal on silicon (LCOS) display, etc.
• CTR cathode ray tube
• LED light-emitting diode
• LCD liquid crystal display
• DLP digital light projection
• VFDs vacuum fl orescent display
• SED surface- conduction electron-emitter display
• Example printing devices may include, for example, toner-based printers, liquid ink-jet printers, solid ink printers, dye-sublimation printers, and inkless printers such as thermal printers and ultraviolet (UV) printers, etc.
• the printing device may print in three dimensions (3-D).
• Output device 175 may receive the object space image information representing the object under observation 125.
• the output device 175 may include, for example, a data storage device, a display device, a printing device, or another computer system.
• the imaging system may further include investigator 170.
• Investigator 170 may be a programmed computer, may receive one or more of the reconstructed image data 261, reconstructed image error data 160, or reconstruction image error data 216, and then apply an algorithm (e.g., pre-programmed routine, artificial intelligence, machine learning, etc.) to extract diagnostic information about the object 125 or to fine-tune control parameters for transmitted s) 115, detector(s) 145, table 130, image reconstructor 165, etc.
• the interface 140 or the output device 175 may not be necessary.
• the investigator 170, the controller 135, and the image reconstructor 165 may reside on the same computer or separate computers.
• Investigator 170 may receive the data 261, 160, or 216 and may be programmed to perform extraction of diagnostic information from the data or to fine tune parameters for processing, for example, at least one of the one or more density domain input pixel, projection directions, projection deformations, projection generation systems process such as focal spot size, etc., or the image reconstructor 165.
• the investigator 170 is also referred to herein as an "investigator computer.”
• Some embodiments of the invention may provide a workstation comprising one or more processors configured to reconstruct an image in a manner similar to image reconstructor 165.
• the workstation may receive input data from at least one of an imaging system, a data storage device, or a computer.
• the input data may be received via a data bus, a cable, a wired network, a wireless network, etc.
• the workstation may further comprise an output device 175 to receive the reconstructed image.
• the output device may be a data storage device, a display device, a printing device, etc.
• Example data storage devices, display devices, and printing devices are as discussed above.
• FIG. 2 shows a systematic chart according to an embodiment.
• FIG. 2A depicts an example raw measurement data input 155, using transformation 200 to compute projection space density data 201 from raw measurements, transformation 205 to obtain transformed projection data 206, and reference transformed prediction 211 using filter 285, together with a projection 276 of the preceding estimated object density 261 for enhanced high efficiency CT with optimized recursions.
• Innovation gain adjustments 220 and 225 form a bridge between FIG. 2 A and FIG 2B.
• FIG. 2B depicts an example forward processing and scaled inversion of the enhanced high efficiency CT with optimized recursions together with the feedback section to obtain the voxel density estimate data 261.
• an invention described herein may include one or more transmitters 115 to transmit an excitation energy 120 into an object under observation 125; one or more detectors 145 to generate projection space data encoding an energy received by the one or more detectors 145 in response to the transmitted excitation energy 120 into the object under observation 125; a controller 135 to control the one or more transmitters 115 to transmit the excitation energy 120 and the one or more detectors 145 to generate the projection space data 155; and an image reconstructor 165 having at least one processor to receive the projection space data 155 and to process the projection space data 155 by, for example, the process depicted in FIG. 2 and explained herein.
• the image reconstructor 165 may: compute projection values for sets of voxels; compute back-projection values for sets of pixels; compute remaining pixels and voxels using other sources of information; and use or build functional relationships among voxels and pixels, such as a priori known high density voxel of a particular object and corresponding expected projection ranges among pixels, or voxels and corresponding missing projection measurements. For example, when some projection values are missing, such as in the missing wedge problem, L. Paavolainen et al. [16], or sparse projections, predictions may be used to replace (expected) missing pixel measurements, using the method of expectation maximization (EM), described in Dempster et al. [17].
• EM expectation maximization
• the EM method represents the case of missing data, equivalent to infinite observation noise.
• Data regarded less reliable (DRLR) in some areas of the image may equivalently be expressed as data with increased observation noise, rather than infinite observation noise.
• DRLR less reliable
• Using all available data, including the DRLR, data in these areas may very efficiently be re-computed/estimated with the EM method, producing re-computed/corrected
• EM data (REMD).
• a weighted combination of the DRLR and the REMD may then be used as input values for the next process iteration, resolving the issue of data weighing dependent on data reliability.
• Relative weighting may, for example, be based on combined noise levels.
• the image reconstructor 165 may re-compute or correct a portion of the at least one unreliable data or missing measurement data using weighted predicted projection data 276.
• An embodiment of the invention may also include, for example, a workstation including one or more processors; and one or more non-transitory data storage devices 180 storing software to be executed by the one or more processors, the software may include software code to implement the process depicted in FIG. 2 and explained herein.
• An embodiment of the invention may also provide a method implemented by one or more processors executing software code stored on one or more data storage devices 180, the method comprising steps to implement the process depicted in FIG. 2 and explained herein.
• the image reconstructor 165 may process the projection space data by the following, or the software may include software code to perform the following, or the method may include steps for the following: transforming the projection space data 155 to obtain projected pixel data 201; nonlinearly transforming the projected pixel data 201 to obtain measured transformed pixel data 206; computing a first pixel innovation result data 216 characterizing a difference between the measured transformed pixel data 206 and predicted transformed pixel data 211; recording the first pixel innovation result data 216 in a data storage device 180; computing a second pixel innovation result data 221, wherein the first pixel innovation result data 216 is re-scaled based on inverses of slopes of two sets of non-linear transformations 205 and 210; computing a third pixel innovation result data 226 based on the second pixel innovation result data 221 from preceding iterations; approximately inverting second pixel innovation result data 221 using a tomographic reconstruction algorithm to obtain preliminary transformed object density update data 236; computing transformed object
• the projection space data 155 is also referred to as "measured pixel data 155," "object imaging data,” “raw measurement data input 155,” and “input data 155.”
• the reconstructed image error data 160 is also referred to as “residual data 160,” “projection residuals 160,” “image data residuals 160,” “innovation residuals 160,” “data 160,” and “residual set 160.”
• the image reconstructor 165 is also referred to herein as “inner feedback loop 165" and "image reconstructor computer.”
• the projected pixel data 201 produced by the cumulative object voxel density along a projection direction using a transformation in box 200, is also referred to herein as "projected integral object density 201,” “projection space density data 201,” “input signal data 201,” “approximate projection data,” “input data 201,” “projections 201,” “data 201,” “projection data pixels set s(i) 201,” and “input projections pixels.”
• the transformation 205 is also be referred to herein as "input transformation functions.”
• the input transformation functions (or input transformation function) 205 may be linear or non-linear.
• the measured transformed pixel data 206 is also referred to herein as “measured transformed pixels 206,” “pixel data 206,” “transformed integral projected object density 206,” “transformed projection data 206,” “innovation data 206,” and “transformed projected object density 206.”
• the transformation 210 is also be referred to herein as "reference transformation functions.”
• the reference transformation functions (or reference transformation function) 210 may be linear or non-linear.
• the predicted transformed pixel data 211 is also referred to herein as “transformed prediction 211,” “predicted data 211,” “predicted transformed pixels 211,” and “transformed prediction data 211.”
• the first pixel innovation result data 216 is also referred to herein as "reconstruction image error data 216,” “image error data 216,” “image data residuals 216,” “innovation residuals 216,” “data 216,” “innovations 216,” “innovation data 216,” “residual data 216,” and “residuals 216.”
• the second pixel innovation result data 221 is also referred to herein as "error data 221,” “innovation data 221,” “innovation residuals 221,” and “data 221.”
• the third pixel innovation result data 226 is also referred to herein as "innovation residuals 226," “data 226,” and “innovation data 226.”
• the preliminary voxel update data 236 is also referred to herein as "preliminary transformed object density update data 236," “preliminary update data 236,” and “preliminary transformed object update data.”
• the transformed object voxel density update data 241 is also referred to herein as "update data 241.”
• the transformed object voxel density update estimate 246 is also referred to herein as “transformed density estimate 246,” “transformed voxel density estimate 246,” and “data 246.”
• the raw voxel density data 251 is also referred to herein as “raw voxel density estimate data 251.”
• the reconstructed image data 261 is also referred to herein as "data 261,” “estimated object density 261,” “voxel density estimate data 261,” “object density data 261,” “reconstructed object 261,” “object data 261,” “object space image,” and “image data 261.”
• MMR-261 Mismatch between a known reference model and the computed object density 261 is referred to as MMR-261.
• the projection 276 is also referred to herein as “predicted projection data 276,” “predicted space data 276,” “predicted projection space pixels 276,” “data 276,” and “predicted feedback measurement pixels 276.”
• the smoothed predicted projection data 286 may also be referred to herein as "data 286.”
• the high processing performance in the image reconstructor 165 may be summarized as the following three step process.
• a feed forward data inversion processing section from projected pixel data 201 to voxel density estimate data 261 that may use data processing linearization at some processing steps allowing efficient use of linear inversion tomographic techniques to capture, ideally, all object related spatial frequency components.
• Nonlinear transformations, allowing computation of small signal voxel gain coefficients, and intermediate transformed object voxel density estimate 246 in order to achieve positive voxel density estimate data 261, are included.
• a feedback loop data projection section computing positive object voxel density estimate data 261 to generate positive predicted projection data 276.
• Projection predictor component H which accounts for the quality of reconstruction by approximating the corresponding object imaging process beginning with the transmitters 115 to data 201, computes an approximation of the object projection process. Optimization with every iteration is highly efficient and based on theoretical properties of the feedback loop rather than on numerical hill- climbing.
• a set of parameters contained in a vector p_, representing uncertainties in an image acquisition system may be chosen and adjusted to correct the effect of these uncertainties. For example,
• uncertainties in the computation of expected scatter or beam hardening may be accounted and adjusted for by computing a corrected vector parameter p.
• the vector parameter p is also referred to herein as "parameter vector p.”
• the adjustments of components of p may, for example, change the base of the logarithmic transformation to account for beam hardening for given object data, applying corrections in the measurement data preprocessing 200, and change the predicted signal intensity data values to account for x-ray scatter using 200 and 275.
• the performance measure may be derived from at least one of a first innovation process 216, projection residuals 160, or expected values of object data of the image reconstructor computer for a fixed set of fixed externally controllable parameter components.
• image data residuals 160 or 216 are sent to the investigator 170.
• the investigator 170 may compute corrected values of parameters in p, using, for example, the highly efficient Levenberg-
• the Levenberg-Marquardt method is typically by a factor k or more times as fast as steepest descent and related approaches, used by others (where k is the number of unknowns to be adjusted). Refined parameters in p are provided back to the image reconstructor
• FIG. 2 may include the functionality of image reconstructor 165.
• measured pixel data 155 of FIG. 1 may be supplied as input as shown in FIG. 2.
• Measured pixel data 155 may be from experimental acquisition 110, simulation 185, or data storage device 180 containing recorded projected pixel data from an earlier experiment or simulation or other data source.
• Measured pixel data 155 may be provided to the image reconstructor 165 directly or remotely via a network, for example, a local area network (LAN), a wide area network (WAN), the Internet etc.
• LAN local area network
• WAN wide area network
• the Internet the global information network
• image reconstructor 165 may process measured pixel data 155 to produce reconstructed object 261.
• measured pixel data 155 is processed to obtain projected pixel data 201.
• the measured pixel data 155 is processed by a transformation that may support the accuracy of subsequent data processing approximations. For example, in x-ray imaging the photon intensity may be transformed by the negative logarithm in order to obtain as a preliminary approximation of an estimate of the corresponding object density projection. It may, for example, contain further refining
• box 200 may use Fourier transformation or principal component methods to convert the measured pixel data 155 into a preliminary approximation of the object density projection space density data, and box 200 may, for example, contain further refining transformations expressed in the parameter vector p that adjust for geometric field distortions and deviations from designed excitation and response measurement conditions.
• the particulars of the transformation process in box 200 may be adapted to approximate any given physical imaging process used to obtain an approximate object density projection density.
• At least a matrix Z (box 220) or a matrix L (box 240) may be determined and used for processing the state data 201 (e.g., the matrix Z (box 220) to obtain transformed pixel innovations 221, or the matrix L (box 240) to obtain voxel updates 241) using the expected state data 276, where the processing depends on the expected state data 276, the function sets of fn (box 205) and f 12 (box 210) to determine the matrix Z (box 220), and the function sets of f 21 (box 290) and g 22 (box 250) to determine the matrix L (box 240).
• the state data 201 e.g., the matrix Z (box 220) to obtain transformed pixel innovations 221, or the matrix L (box 240) to obtain voxel updates 241
• the processing depends on the expected state data 276, the function sets of fn (box 205) and f 12 (box 210) to determine the matrix Z (box 220), and the function sets of f 21 (box 290)
• projected pixel data 201 containing information encoding the internal structure of an object may be non-linearly (or linearly) transformed into pixel data 206.
• f 11 may be computed in box 205 as an approximately variance stabilizing function of the input data 201.
• f 11 in box 205 may be approximately the Poisson variance stabilizing square-root function of the input signal data 201 obtained from PET measurements.
• f 11 may be a linear function or a set of linear functions, for example, when the measurement noise is constant.
• predicted projection data 276 may be transformed into one or more predicted transformed pixels 211 using nonlinear (or linear) reference transformation f 12 pixel- by-pixel.
• pixel may refer to the location of the pixel or the value of the pixel.
• voxel may refer to the location of the voxel or the value of the voxel.
• the data is transformed from the density domain to the transformed domain.
• the nonlinear (or linear) reference transformation f 12 may use one or more predicted projection space pixels 276 (or smoothed predicted projection data 286 when smoothed by a filter (such as a low pass filter) in box 285) from the density domain to compute one or more predicted transformed pixels 211 in the transformed domain.
• a filter such as a low pass filter
• The, for example, two- dimensional smoothing characteristics, in box 285 on a projection 276 of an object are an approximation to the effect of the projection when smoothing the same object data with a filter 270 in, for example, three dimensions. This concept also applies to higher dimensional reconstructions. In other words, the linear operations of projection and smoothing are interchangeable.
• the predicted projection space pixels 276 may be obtained from the predictor H (in box 275), using the prior voxel density estimate data 266.
• the nonlinear (or linear) reference transformation f 12 will tend to be pixel-by-pixel an approximation of f 11 .
• f 11 and f 12 may differ, for example, when elements of f 11 require adjustments for artifacts in data 201, to allow robust estimation processing.
• f 12 may be a linear function or a set of linear functions.
• innovation data 216 may be calculated in the transformed domain.
• one or more measured transformed pixels 206 may be compared with predicted transformed pixels 211, the non-linearly transformed data from smoothed predicted projection data 286 to produce a difference corresponding to innovation data 216.
• the difference between corresponding measured transformed pixels 206 and predicted transformed pixels 211 may be calculated to produce the first pixel innovation result data 216.
• the first pixel innovation result data 216 may be used for residual analysis and for optimizing systems parameters to ultimately yield the best possible object density reconstruction.
• each of the corresponding pixel locations weighted combinations z of the inverse slopes of the nonlinear input transformation fn from box 205 and the inverse slope of the nonlinear reference transformation f 12 from box 210 of the associated density domain pixel values (in the neighborhood of the corresponding pixel locations) may be individually computed.
• the slopes corresponding to particular density domain pixel values may include their spatial and their sequentially neighboring values, and may be subject to reevaluation using robust functions.
• first pixel innovation result data 216 may be scaled by a pixel -by-pixel matching matrix Z to generate error data 221.
• the pixel -by-pixel matching matrix Z may be computed from the current projection space data and characteristics of their transformation.
• the scaling matching matrix Z in box 220 may be calculated by organizing a set of z's, such that each z corresponds to an associated pixel location. The calculation may use inverse slopes of the nonlinear input transformation functions 205 and inverse slopes of the nonlinear reference transformation functions 210 for corresponding input and reference pixel values.
• the coefficients of matrix Z may be representations of the combined inverse slopes of the pixel transformation function fn at pixel values of one or more measured transformed pixels 206 and fo at pixel values of predicted transformed pixels 211.
• the combination of the inverse slopes z uses a process that is likely to produce near optimal values z-opt in the noiseless case and combines them to be likely to produce a close approximation to z-opt in the noisy case.
• a second pixel innovation result 221 may be computed using the pixel -by-pixel product of the first pixel innovation result data 216 and the corresponding elements of the pixel innovation scaling matrix Z.
• the pixel-by-pixel innovation scaling matrix Z may be calculated using for each pixel an inverse slope of the corresponding nonlinear input pixel density transformation function 205 and an inverse slope of the corresponding nonlinear reference pixel density transformation function 210.
• the nonlinear input pixel density transformation function 205 may be range limited or robust (e.g., constrained values).
• the nonlinear reference pixel density transformation function 210 may be range limited or robust (e.g., constrained values).
• a matrix gain refinement equivalent to the technique of U.S. Patent No. 8,660,330 based on the sequential properties of the pixel and neighboring first pixel innovation result data 216 may be optionally integrated.
• a second in-series innovation gain matrix G equivalent to the technique of U.S. Patent No. 8,660,330 based on the sequential properties of the pixel and neighboring pixel innovation results may transform the second pixel innovation result data 221 into a third pixel innovation result data 226.
• elements of the pixel-wise correcting gain matrix G may be associated with a corresponding pixel in the second pixel innovation result data 221 and the third pixel innovation result data 226.
• the second pixel innovation result data 221 may be used to compute the third pixel innovation result data 226 using, for example, the gain evaluation technique based on the sequential properties of the first pixel innovation result data 216 or second pixel innovation result data 221, corresponding to the gain evaluation technique of U.S. Patent No. 8,660,330.
• the raw corrective gains may be determined corresponding to the gain evaluation technique of U.S. Patent No. 8,660,330, using the sequential properties in any of the innovation processes.
• the pixel-wise correcting gain matrix 225 is computed by regression on one or more spatially neighboring coefficients of innovation gains or one or more sequential coefficients of innovation gains from preceding iterations.
• the correcting gain matrix 225 represents measurement and model defects.
• the raw correcting coefficients prior to regression are computed in an innovation processor 230 using at least one of the first pixel innovation result data 216 or the second pixel innovation result data 221, to identify innovation patterns A or B, and act with corresponding raw correcting increases or raw correcting decreases of coefficients.
• the third pixel innovation result data 226 may be provided to a tomographic reconstruction algorithm to obtain preliminary voxel update data 236.
• the tomographic reconstruction algorithm may be, for example, a linear reconstruction algorithm
• the tomographic reconstruction inversion algorithm using the DC/average projection value of Zeng et al. [1] may be used.
• the tomographic reconstruction algorithm may be a back-projector.
• the relations between the necessary projection DC/average values and other frequency components, characterized by a set of systems parameters, may depend on the amount and orientation of available projection data. For example, for two (noiseless) orthogonal projections, reconstruction with the invention may not require filtering of projections for reconstruction.
• the preliminary transformed object density update data 236 may be rescaled with the elements of a voxel -by- voxel update scaling matrix L, where each transformed object voxel density update data 241 may be associated with an element of the matrix L.
• the vox el -by-vox el update scaling matrix L may be computed from the transformed object data and the characteristics of a positive constraining transformed object voxel data transformation.
• Values of the elements of the matrix L may be computed from the inverse slope of the elements of the voxel-by-voxel back-transformation using voxel-by-voxel elements of g 22 in box 250 of the corresponding transformed voxel density estimate 246.
• the transformed density estimate 246 may be calculated by adding the transformed object voxel density update data 241 to the corresponding transformed preceding voxel data estimate 291.
• the transformed density estimate 246 may be stored in data storage device 180 for later use by the investigator 170.
• the transformed density estimate 246 may be back-transformed into the density domain to obtain the raw voxel density estimate data 251.
• the data is transformed from the transformed domain to the density domain.
• Transformed density estimate 246 may be voxel-by-voxel non-linearly transformed using the matrix elements g 22 to produce raw voxel density estimate data 251.
• the matrix elements g 22 may be function types as an approximate inverse of the, one or more, function types of the matrix elements of fn or f 12 for major sections of the matrix g 22 , and the matrix elements of g 22 may be used to compute output values satisfying external inputs in at least one region of the data range of the transformed density estimate 246.
• the functions of box 250 are also referred to herein as a "positive constraining functions.”
• the raw voxel density estimate data 251 may be filtered to obtain the preliminary voxel density estimate data 256.
• Box 255 may use a low pass filter to smooth the raw voxel density estimate data 251.
• Box 255 may filter the raw voxel density estimate data 251 to output the preliminary voxel density estimate data 256 as positive.
• the low pass filters in boxes 270 and 285 may share the property of keeping their outputs positive, and may replace the function of box 255.
• the preliminary voxel density estimate data 256 of the object may be post-processed by P to operating on a single voxel at initial iteration and a multiplicity of voxels subsequently, refining density values using suitable a priori information about the properties of voxel density estimate data 261.
• This technique may start with a single voxel, but does not require a plurality of voxels at the outset.
• Post-processing the preliminary voxel density estimate data 256 may create one or more voxels based on one or more of (a) equal or increased grid resolution or (b) equal or increased density value resolution.
• voxel density estimate data 261 may be recorded, delayed, and transmitted as prior voxel density estimate data 266 as inputs to the filter of box 270 and then to the transformation f 21 of box 290, and also as inputs to the predicting projection processor H of box 275.
• the prior voxel density estimate data 266 is filtered to obtain smoothed prior voxel density estimate data 271.
• the filter of box 270 may be a low pass filter.
• the data is transformed from the density domain to the transformed domain.
• the nonlinear reference transformation family f 21 may also be referred to herein as a feedback function 290 and as a set of feedback functions in box 290.
• the prior voxel density estimate data 266 may be processed to obtain predicted projection data 276, typically using a representation of a single sparse projection matrix H.
• the matrix H may be designed to duplicate the imaging of the original object density values as object projection density values. Uncertainties in the imaging process in FIG. 1, however, may result in uncertainties in the matrix H and may be expressed in a vector parameter P for adjustment when more is learned about the imaging process.
• the matrix H may express scatter for a particular voxel using patterned off-voxel-projection components, parametrized in the vector p. For practical processing, however, explicit use of the matrix H may be avoided.
• box 275 functionality may be related to those discussed in, for example, Zhang et al. [2] and Long et al. [3],
• the projection process in box 275 summarized in the matrix H may, for example, represent at least one system parameter or object parameter.
• the system parameter may be, for example, one or more of the focal spot geometry, focal spot beam exit intensity and hardness exit profile characterization, changing tube supply voltage, beam dependent x-ray detector characteristics, or x-ray scattering.
• the object parameter may represent, for example, object movement.
• a tomographic projection algorithm may be used with a vector parameter P in H(p). After convergence of the processing loops shown in FIG. 2, and without changing p ⁇ the parameter p may be adjusted to search for improved solutions to performance criteria.
• the performance criteria may, for example, be based on the final set of residuals of the innovation process 216, the residual set 160, or a mismatch of values between the computed and reference object values.
• the nonlinear transformation f 12 may be applied to the predicted projection data 276, obtaining the transformed prediction data 211.
• Vector parameter p uses data 216, 221, 226, or 160, for example, for their weighted sum-of-squares minimization, to optimally fit systems parameters to changing projection conditions. Projection conditions are expressed, for example, in parameters of the projection operation in H in box 275, post processing method P summarized in box 260, and parameters in functions fn in box 205, functions f 12 in box 210, functions f 21 in box 270, and functions g 22 in box 250.
• Data 246 and 261 may be modified to satisfy external systems variables and objectives. For example, a resetting of some object density variables contained in data 246 or 261 to fixed a priori known values may support more rapid convergence of remaining object density variables and systems parameters.
• the predicted projection data 276 may be computed based on a set of one or more object space voxels, wherein the set of object space voxels may cover a plurality of resolution grids (e.g., including a single grid point) of varying sizes when the set of the object space voxels are projected onto projection space.
• the set of object space voxels may cover a plurality of resolution grids (e.g., including a single grid point) of varying sizes when the set of the object space voxels are projected onto projection space.
• the residual data 160 may be produced from the projected pixel data 201 and the predicted projection data 276.
• the residual data 160 may be transferred to the investigator 170 and the output device 175.
• a correcting gain matrix may be computed at each innovation using a spatial and sequential weighting function weighing one or more raw, neighboring and preceding, innovation gains computed with a method that uses, for example, at least some of the features of the loop gain determination used in U.S. Patent No. 8,660,330. Implementing features of U.S.
• Patent No. 8,660,330 to the present feedback loop of FIG. 2 increases raw gain coefficients in box 230 when innovations 216 or 221 are consecutively statistically significant of equal sign, creating pattern A, and decreases raw gain coefficients in box 230 when innovations are consecutively statistically significant of alternating sign, creating pattern B.
• Different versions of this approach are possible. A more primitive process version would be, starting from
• the invention may compute a correcting gain matrix coefficient using a weighting function that may be computed by regression on one or more spatially neighboring innovation matrix coefficients or one or more sequential innovation matrix coefficients from preceding iterations to account for measurement and model defects.
• the one or more sequential innovation matrix coefficients may also be referred to as one or more sequential coefficients of innovation gains.
• An object of interest may be refined iteratively by the investigator in box 170 controlling p_, using the weighted and converged innovation residuals of 160, 216, 221, or 226 created by final iterations for each set of p_ within the loop of FIG. 2.
• Perturbations of the parameter p_ computed in the outer loop, where each set of weighted innovation residuals is taken from sufficiently well-converged iterations within the inner feedback loop 165, shown in FIG 1, and summarized in FIGs. 2 A and 2B, may be used for optimization.
• the externally controllable components of the parameter p_ may gradually change from an initial estimate to that forming the best match to external objectives or minimal cost in terms of innovation residuals 160, 216, 221, 226, or, for example, deviations from a known systems reference model, such as in simulations, or when evaluating convergence of an object density estimate to the density of a reference object.
• the investigator computer [170] externally to the image reconstructor computer [165] may be operable to gradually change controllable components of a parameter vector from an initial setting to those producing the lowest cost relative to an external objective, while using properties of at least one of the projection residuals or object residuals characterizing the difference between a known reference object and the corresponding computed object density.
• Projection residuals may be processed with robust estimation influence functions that have been vetted through simulation of object and recording challenges using a priori parameters or parameters that represent and/or are derived using their neighborhoods.
• the parameters may be, for example, functions of one or more of local input measurements, systems properties, prediction values, or a priori expected object and measurement properties.
• Projection residuals may be processed with at least one of smoothing the residuals over select ranges, scaling the residuals over select ranges, or weighing by an influence function over select ranges.
• Innovation data 216, 221, 226, preliminary update data 236, and update data 241 may be processed using influence functions.
• Measured pixel data 155 may be processed using influence functions.
• mismatch residual MMR-261 may be processed using influence functions.
• Embodiments of the process in FIGS. 2 A and 2B may include, for example, extended high efficiency computed tomography with optimized recursions (eHECTOR) that uses the ability of a linear reconstruction algorithm (LRA) (e.g., filtered back-projection (FBP) or other qualified, preferably linear, reconstruction algorithm) to address eigenvalues of the tomographic inversion problem efficiently.
• LRA linear reconstruction algorithm
• FBP filtered back-projection
• the LRA may be embedded in a non-linear structure that may be linearized using pixel-by-pixel small-signal data gains in 220, 225, and
• FBP may be used here recursively on a linearized model of a transformed estimation problem.
• This linearized loop of FIG. 2 represents a highly contractive mapping inducing geometric shrinking of estimation errors. This shrinking will terminate when residual errors, such as measurement noise, or model errors induce inconsistencies that may not be resolved by further iterations. Inconsistencies may be seen in projection residuals and are a driving force for re-estimation of the vector parameter p.
• the LRA is a suboptimal approximation for the optimal filter gain.
• the sequential processes of iterations in box 165 neglecting the low-pass data filter in box 255, the increase of grid resolution with small innovations, and the use of an approximate, sub-optimal, not explicitly computed LRA gain, may be compared with an optimal linear Kalman filter with fixed dimensionality of its state variables, explicitly computed gain K, and well specified prior statistics with finite variance.
• Table 1 depicts a discrete, optimal, Kalman filter algorithm, edited, as shown in Sage et al, p. 268 [4].
• Uncorrelated observation noise V v ( i) is assumed here despite using the same fixed set of measurement components, rather than a new set of independent measurements-noise contributions. Sequential correlation is neglected, for example, in view of smoothing operations, similar increases of grid resolution, and redundant measurements, because noise data may then be mutually inconsistent. As such, sequential correlation (almost) cannot be expressed in the object density (except for the damaging image noise) and the state update may therefore use the same computational step as in the presence of white noise. This model may be justified to the extent that reconstruction of a fixed set of noise values without the presence of an object 125 will produce negligible voxel densities, when compared to an eventual presence of the object 125.
• boxes 255, 260, 270, and 285 may suppress high-frequency components in voxel data values, leaving correspondingly large residual projection values due to, for example, "object projection inconsistency.”
• Vw( 0) represents the initial object density variance contained in the state variance Vx(i)
• H( i) the observation matrix for the object density state vector x( i) and [1.5] v( i) is the observation / measurement noise. [1.6]
• K( i + l) K( i) [1.8] during iterations with fixed grid resolution and while determining the operating points for data processing linearization in the corresponding eHECTOR.
• FIG. 2 also shows a special case of using, for example, logarithmic functions fn (box 205), logarithmic functions fi 2 (box 210), logarithmic functions f 2 i (box 290), exponential functions g 22 (box 250), innovation scaling matrix Z (box 220), and update scaling matrix L (box 240).
• both matrices Z and L respectively, may be represented by diagonally dominant matrices, each diagonal element corresponding to data elements such as a single pixel or a single voxel in the data stream.
• the elements of matrix Z(i) may represent pixel-by-pixel a weighted average of a first component of small-signal inverse slopes of function fn(s(i)) (box 205) for projection data pixels set s(i) 201, and a complimentary corresponding weighted second component of the small-signal inverse slopes of function fi 2 (pf(i)) (box 210), where pf(i) may be the set of values of predicted feedback measurement pixels 276.
• the elements of matrix L(i) may represent voxel-by-voxel the small-signal inverse slopes of function g 22 (dtr(i)) or its weighted average with the small- signal slope of function f 2 i (266).
• the function fn (box 205) and the function fi 2 (box 210) may form a pair of identical sets, while the function f 2 i (box 290) and function g 22 (box 250) may form a pair of sets of mutually inverse functions.
• the functions for functions fn (box 205), fo (box 210), and f 2 i (box 290) may be, for example, the logarithmic function (for example, for very small projection count) or the square-root function (for, for example, photon count signal variance stabilization of measured positron emission tomography (PET) or single photon emission tomography (SPECT).
• the functions g 22 may be, for example, the exponential or portions of a quadratic function.
• the function fn (box 205) and the function f 12 (box 210) may be, for example, a Poisson variance stabilizing square-root function of the one or more input projection pixels.
• the generalization adds matrices Z (box 220) and L (box 240), related to the non-linear function choice for fn (box 205) and f 12 (box 210), and f 21 (box 270) and g 22 (box 250), shown in FIG. 2.
• Table 1 yields for the static linear model shown in FIG. 2:
• a new Kalman gain K' may be specified and associated with the small signal loop gain-matrices Z (box 220) and L (box 240).
• the small signal loop gain-matrices Z and L may be based on the slopes of the nonlinear transforms.
• gain matrix Z (box 220) may be used to compensate for the slopes of the functions fn and fi2 associated with the operating point of the input transform in box 205 (where the projected pixel data 201 comes in (e.g., logarithmically transformed) and the prediction data 211 from box 210.
• gain matrix L (box 240) may be associated with the output transformation g 22 (e.g. exponential) and its inverse f 21 .
• the Kalman K matrix (box 235) may be, for example, approximated by a LRA such as the filter back-projection (FBP) or the back-projection filter (BPF) approximation (see, e.g., Zeng et al. [1]), and H may be the forward projection of the estimated object, including phenomena such as scatter.
• LRA such as the filter back-projection (FBP) or the back-projection filter (BPF) approximation (see, e.g., Zeng et al. [1])
• H may be the forward projection of the estimated object, including phenomena such as scatter.
• the optional gain matrix G (box 225) may support corrections to K (box 235) because K is usually approximated, and G (box 225) aids to compensate effects such as the entropy increasing, stabilizing, and constraining low-pass filter, a priori constraints, and other interventions in the feedback loop.
• the low-pass filter 255 and a priori object knowledge 260 may be set up to support stability of eHECTOR, especially in sparse projection data object density computation.
• the small signal analysis may start with given input data 201 for the degenerate case of a single grid point for all matrices and adjusting, for example, scalar gain (matrix) G to provide unity (DC) loop gain.
• the initial loop gain may, for example, be adjusted using the matrix G, based on the observed innovation sequence data similar to U.S. Patent No. 8,660,330, or other criteria.
• the use of the smoothing operation (entropy increase operation) in box 255 may be neglected or adjusted to have a minimal effect.
• the non-linear filtering problem in FIG. 2 may be analyzed using the Kalman filter equations shown in, for example, Table 1.
• the small signal diagonal gain matrix Z (box 220) may have the following inverse slope elements (assuming
• Zd represents the diagonal non-zero Z-matrix element with index d
• ⁇ is the differential operator.
• the slope zd may also be a weighted combination derived from the function fi2 and the function fn depending on noise properties of the data 201 and 286.
• the use of coefficients zd may be subject to constraints by influence functions which are discussed further below.
• the approximated small signal perturbation may then be passed through the LRA (box 235), approximating the Kalman filter gain, computing the preliminary voxel updates.
• the small signal output of the filter (data that corresponds to x( i) in Table 1) may be followed by a small-signal transformation matrix L (box 240) compensating the subsequent nonlinear gain of the transformed object density estimate dtr ( i).
• the diagonal matrix L (box 240, relative to voxel density v q ) has elements: replacing index qq ⁇ q to indicate the diagonal elements.
• the use of coefficients l q may be subject to constraints by influence functions which are discussed further below.
• a new small-signal data perturbation projection matrix H' is obtained, defined as:
• Recursively passing through the feedback loop yields for the small-signal approximation a geometric object density contraction factor limited by the deviation C x of the approximation of the LRA to the optimal Kalman filter gain.
• a geometric object density contraction factor limited by the deviation C x of the approximation of the LRA to the optimal Kalman filter gain.
• the a posteriori variance algorithm for the linear Kalman filter in Table 1 showing the reduction of posterior variance as a function of the match between I and KH.
• inconsistent measurement components lead to weighted least-squares estimation without further contractive mapping.
• the structure of the residuals at convergence aid in diagnosing systems performance. For example, white noise projection residuals and negligible object features, if any, in the projection residuals may indicate convergence and quantify accuracy of the model obtained in reconstruction.
• Judgment on C x ( i) expressing the efficacy of K( i) may be based on:
• K may be augmented with an innovation gain adjuster matrix G (box 225).
• the gain adjuster G (box 225) modestly modifies gain values in the range of the operating points for the slopes of fn of the input data (box 205) and the slopes fi 2 of the prediction data of the observations (box 210).
• the gain adjuster matrix G may be influenced by factors such as beam hardening, resulting from consistent under estimation or over estimation of predicted projection density.
• the value of using the matrix L is to compensate the small- signal loop gains associated with the functions g 22 supporting an overall small signal loop gain from 216 to 211 of the identity I (equation 2.9), their difference forming in box 215, implying a small difference aside of noise, and consequently rapid convergence of the object density during iterative computations.
• the function fn (box 205), the function fi 2 (box 210), the function f 2 i (box 290), and the function g 22 (box 250) may be used for matrix Z in box 220 and matrix L in box 240
• some of the small signal matrix elements may be set to zero or be left out in the presence of some a priori information.
• the object density may be known a priori at some points, such as for implants or the density surrounding the patient.
• Elements associated with estimating these object densities may be left out, and corresponding object density data replaced with the known data values.
• parameters that are not known or observed by other means expressing, for example, uncertainty in beam hardening due to object properties such as implants, beam spectral properties, variable emitted beam hardness, object movement, scatter, respiration and cardiac movement, and other systems components outside of the image reconstructor 165, may be represented and adjusted inside the image reconstructor 165, for example, in the projection matrix H (box 275), shown in FIG. 2.
• object properties such as implants, beam spectral properties, variable emitted beam hardness, object movement, scatter, respiration and cardiac movement, and other systems components outside of the image reconstructor 165
• Such parameters may also adjust the functions in boxes 200, 205, 210, 250, 255, 260, and 290.
• These parameters may be estimated jointly or separately, for example, using the Levenb erg-Mar quardt (LM) approach.
• LM Levenb erg-Mar quardt
• the LM, or equivalent, approach computes from the residual data, which in turn were computed from sets of changed parameters inside the image reconstructor 165 when their corresponding object estimates have converged, adjustments to the parameters with uncertainty continue until they have converged.
• FIG. 3 and FIG. 4 An example of such LM-based reconstruction is shown in FIG. 3 and FIG. 4, using the LM processing for the estimation of a beating model of a coronary tree, the data being collected with a clinical bi-plane C-arm system.
• the model parameter computation used over 50 parameters to re-estimate the alignment of projections resulting, for example, from static and variable magnetic field image intensifier distortion, and C-arm position errors.
• the initialization of the iterations in image reconstructor 165 may start at a lower grid resolution, for example from a single voxel and a single pixel for each projection, rather than a multiplicity of voxels.
• the initial voxel density follows, for example, from the average of the projected densities.
• the grid resolution may be increased to the desired level. Changes of resolution may not be at equal iteration intervals or with equal resolution scale-steps.
• the number and sequence of incorporating components of model parameters p_ as grid resolution is increased may be predetermined, based on reduction of errors, or information criteria as in, for example, Schwarz [6] and Akaike [7].
• Acceptance of parameter coefficients for example, may be based on measures of data 206, or 160, and the degree of improvements incorporating components of p_ for reconstruction, and may be used for determining the importance of parameters and their selection.
• the major gain adjustments of the loop gain are determined by box 220 and box 240.
• Remaining gain adjustments may be implemented using the gain in matrix G(i) (box 225) that may be derived from the patterns of sequences of preceding innovations 216.
• Significant oscillations at fixed grid resolution may indicate excessive gain, persistent and significant static innovation values insufficient gain, and random patterns appropriate gain (see, e.g., U.S. Patent No. 8,660,330).
• Nonlinear iterative tomographic reconstruction has been pursued intensively by many and for many years.
• Elbakri et al. [8] and [9] describe approaches to account for the polyenergetic X-ray source spectrum, energy dependent attenuation, and non-overlapping materials. Their objective is to compute an unknown object density using the known energy- dependent mass attenuation coefficient.
• Elbakri's approach is to formulate a penalized- likelihood function for this polyenergetic model and using an ordered subset iterative algorithm to estimate the unknown object density for each voxel.
• Staub et al. [11] demonstrated the ability to use and track the motion of the patient anatomy on a voxel by voxel scale, using a small number of eigenvector representations for object reconstruction.
• Staub et al. [11] found the Nelder-Mead simplex algorithm the most robust approach to estimate the small number of parameters representing object movement.
• Reconstruction parameters may, for example, use the residual sum of squares or other criteria of computed robust estimates of residual data using influence functions, for example of data 160, data 216, and, for example in the case of simulations, mismatch residuals (MMR-261), to determine reconstruction performance.
• influence functions for example of data 160, data 216, and, for example in the case of simulations, mismatch residuals (MMR-261)
• MMR-261 mismatch residuals
• optimal object reconstruction is performed for any set of parameters, equipment and object. From the parameter vector p_ a globally optimal set of parameter data p_ is computed, that minimizes reconstruction errors or optimizes reconstruction performance.
• the invention described herein combines highly contractive parameter estimation techniques (equation 2.9).
• the effectiveness of the invention here may be derived from, for example:
• Parameters may represent, for example,
• FIG. 6 Exemplary results of this approach are shown in FIG. 6.
• the innovative eHECTOR demonstrates higher image quality and more than 10-fold higher speed when compared with CERN's August 30, 2017 release of SART (Biguri et al. [21]).
• FIG. 6A eHECTOR shows significantly more precise object feature reconstruction and smaller residual error than SART.
• FIG. 6B the preliminary version of eHECTOR demonstrates useful object reconstruction results from sparse projection data sets, e.g. five projections, with those obtained by CERN's SART (32 iterations).
• 6A and 6B show eHECTOR and SART implemented at the Basser Laboratory of the Eunice Kennedy Microver National Institute of Child Health and Human Development (NICHD), National Institutes of Health (NIH). This is not an endorsement of the eHECTOR process by the NIH.
• model parameter fit optimization with respect to a subset of the parameter vector p, controlled, for example, by investigator 170 in FIG. 1, minimizes these residuals representing the error between measurement and prediction.
• weighting When using multiple measurement methods, such as dual energy scanning, their weighting has to be specified relative to the cost functions of the estimation problem.
• This numerical optimization may be used at any given analytics-based converged iteration during the iterative reconstruction algorithm.
• a preferred method for numerically estimating p subsequently, is the Levenberg-Marquardt procedure
• LMP Long ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇ ⁇
• NLR near-linear regions
• termination of any of the iterations, including at the highest resolution level, may be predetermined, limited, or derived, for example, from statistical measures and other cost criteria reaching thresholds.
• variable object density beam hardening absorption coefficients can model functions of local object density and the effect of the integration path-length on projections, as desired for Elbakri's et al. [8] or
• Humphries' estimation objective may also be used to estimate parameters characterizing the imaging system, such as, for example, focal spot size, misalignment of projection directions, and advance of the patient table during the image acquisition.
• this method may be preferable to the use of fixed prior reconstruction parameters allowing instead self-calibration of the system, "life" systems performance validation, and reduced maintenance service requirements.
• SNR signal-to-noise ratio
• MIR Model Based Iterative Reconstruction
• PWLS penalized weighted least squares
• a variance stabilizing transformation Sv in box 205 simplifies estimation and systems performance evaluation. For some measurement distributions computational transformations can be used based on known characteristics; for others, transformations may be evaluated numerically, rather than theoretically.
• stabilizing transformations Sv for box 205 are variations of the square-root transformation.
• An alternative approach to reduce the influence of noise Vv on object density variability is to combine projection measurements locally, e.g. by local smoothing, comparable to a reduction of the dimensionality of the estimation.
• Influence functions can reduce the effect of outliers or departures of measurements from model assumptions and allow improved object reconstruction with little computation and little quality loss when compared to data without defects.
• Influence functions may be of use in any of the connections between boxes in FIG. 2. Most likely, however, influence functions may be located between boxes 200 and 205 or between boxes 215 and 220 and boxes 235 and 240 or between boxes 240 and 245.
• the influence function can be designed to account for data with variable SNR.
• the influence functions restrict maximal amplitudes of signals and may be redescending beyond critical input values.
• Bouman et al. [13] study the properties of influence functions that preserve edges, and show their usefulness in maximum a posteriori (MAP) and log-likelihood estimation. Their objective, however, is to retain certain features in random fields. This is in contrast to the invention, where consistent measurement structural information may be accumulated in the estimation of the object density and randomness / entropy may be expressed in residuals.
• a first influence function between boxes 215 and 235 may, for example, be applied to data 216 or contained in box 220 by modifying the coefficients z, or modify the gain coefficients in box 225, operating on the residual data following the comparison of measurement and reference signals in box 215, but prior to the inversion that generates the object density update data signal in box 235.
• the prime candidate subject to the influence function is data 216, the residuals between transformed measurements and transformed prediction.
• the choice of parameters specifying this influence function depends on the values of its corresponding joint neighborhoods in pixel data 206 and 211, the distribution of function values within the pixel neighborhoods of data 216, and its change relative to previous iterations. Location and shape of any influence function within FIG. 2 may tend to be problem-specific and result from the field of the application.
• FIG. 5 shows an example of a shape of an influence functions for use on data 216 respectively in box 220.
• the influence function depicts a processing dependence on data indicating the reliability of measurements, such as associated with data 201, 206, or 211.
• the influence function may make subsequent data processing robust against, for example, rare large erratic measurement in an environment of noisy measurements.
• Asymptotic slopes may be 0.00.
• the influence function can be part of the innovation gain adjustment in box 220.
• reconstruction of dense masses may lead to low SNR of their estimates due to corresponding low photon counts and relatively large variations of the associated projected pixel data 201.
• the simplicity of using influence functions, as shown in FIG. 5, for example, in combination with variance stabilization is computationally preferable to alternative approaches using the stabilizing penalty function S(x) with the numerically demanding minimization of the model likelihood D(y; x), when
• Exemplary embodiments of the invention may be provided as software code stored on a non- transitory data storage device, such as, for example, CD-ROM, DVD, BLU-RAY, magneto- optical (MO) disk, hard disk, floppy disk, zip-disk, flash-drive, etc.
• the stored software code may be readable and executable by a computer having one or more processors, one or more non- transitory memory devices, such as, for example, random-access memory (RAM) devices, dynamic RAM (DRAM) devices, flash memory devices, and static RAM (SRAM) devices, etc., to perform the exemplary techniques discussed above with respect to, for example, FIG. 1 and FIG 2.
• RAM random-access memory
• DRAM dynamic RAM
• SRAM static RAM
• Exemplary embodiments of the invention may provide one or more program storage and execution devices, for example, application specific integrated circuits (ASICs), field programmable gate arrays (FPGAs), complex programmable logic devices (CPLDs), application specific instruction-set processors (ASIPs), etc. for storing and executing the exemplary techniques as discussed above with respect to FIG. 1 and FIG. 2.
• ASICs application specific integrated circuits
• FPGAs field programmable gate arrays
• CPLDs complex programmable logic devices
• ASIPs application specific instruction-set processors

Priority Applications (2)

Applications Claiming Priority (4)

Publications (1)

Family

ID=62146811

Family Applications (1)

Country Status (4)

Cited By (1)

Families Citing this family (7)

Citations (3)

Family Cites Families (15)

• 2017
  • 2017-11-17 JP JP2019525981A patent/JP7280184B2/ja active Active
  • 2017-11-17 EP EP17872842.4A patent/EP3526733A4/en active Pending
  • 2017-11-17 US US15/816,745 patent/US10607378B2/en active Active
  • 2017-11-17 WO PCT/US2017/062369 patent/WO2018094250A1/en unknown
• 2019
  • 2019-11-25 US US16/694,601 patent/US10789743B2/en active Active

Patent Citations (4)

Non-Patent Citations (6)

Also Published As

Similar Documents

Legal Events