Radiation scanning of objects and, more particularly, radiation scanning of objects, such as cargo containers, by multiple sources, to reconstruct three-dimensional images.
Radiation scanning systems used in cargo imaging technology, including land, sea, or air cargo, have advanced through multiple generations. First generation radiation scanning systems include gamma-ray sources, which are used in the mobile VACIS® gamma ray imaging system, available from Science Applications International Corporation, San Diego, Calif. (“SAIC”). The mobile VACIS® includes a Co-60 or Cs-137 isotope source and an array of individual NaI detectors operating in photon counting mode, mounted to a truck. Due to safety concerns and finite detector response time, source strength is around 1 Ci. Imaging resolution, penetration and contrast sensitivity are poorer than necessary. Typical penetration is 100 mm to 150 mm steel with ˜15 mm display resolution (pixel size).
Second generation radiation scanning systems use single energy X-ray sources, such as the Varian M-series Linatrons, available from Varian Medical Systems, Inc., Palo Alto, Calif. The use of intense pulsed X-ray source and integration detectors greatly improved radiographic performance, such as resolution, penetration and contrast sensitivity. Steel penetration is over 400 mm with 3 mm resolution with an M6/Mi6 source. This allows an inspection to reveal most (two dimensional) details inside a cargo container.
Third generation scanning systems use dual-energy X-ray transmission radiography with material discrimination. In addition to superior radiographic performance, this generation of X-ray scanning systems provide pixel-by-pixel information of material classes, such as organic, inorganic, metallic and very high-Z metals. Algorithms may also be applied to smooth out material regions based on the fact that material types is not likely to change from pixel to pixel.
Fourth generation radiation scanning systems include dual-energy X-ray transmission radiography with stacked detectors and data fusion algorithms for enhanced material classification. While stacked detectors by themselves may not work well at MeV energies, they pick up material signatures at small beam path length, which complements information from dual source energies. Coupled with advanced algorithms, fourth generation radiation scanning systems can provide better material classification than the prior generations discussed above, and is useful with small beam path length as well as larger large objects.
The linear detector array 18 is electrically coupled to an image processor 20, which is coupled to a display 22. The image processor 20 comprises analog-to-digital conversion and digital processing components, as is known in the art. A processing device 24, such as a computer, for example, is electrically coupled to and controls the operation of one or more of the electron accelerator 12, the linear detector array 18, a conveyor system (not shown), the image processor 20, and the display 22. One or more memory devices 26 to store a reconstruction algorithm, detected data, resulting images, etc., is also provided. Connections between the processing device 24 and of all the components are not shown, to simplify
Conventional CT reconstruction typically requires hundreds or thousands of views. Filtered backprojection (“FBP”), or, more broadly, direct or analytic algorithms are commonly used to reconstruct CT images. Such algorithms are derived as analytic solutions to some idealized version of the actual problem to be solved. Two problem aspects that are commonly idealized are the geometry and the physics. Geometrically, most analytic algorithms assume that data are continuously sampled along an ideal arc or line. In terms of physics, most also assume an ideal X-ray imaging mechanism, such as a noise-free monochromatic infinitesimal-width pencil-beam from a point source with no scatter using perfect noise-free linear electronics, for example. Both of these assumptions are problematic.
Known fourth generation systems (dual-energy source with stacked detectors and algorithm) only provide projection images. Objects extending through an X-ray beam path are “squeezed” into a planar view. Except in a few simple cases, the 2-D image (with material information) may not contain sufficient information to resolve objects that overlap in a beam path direction. Furthermore, material information (atomic number or Z-value) is calculated for all materials along each beam path as a whole. This dilutes the material signature. A large piece of ordinary material will reduce the signature of a heavy metal along the same beam path.
In full scale X-ray tomography, projections, which are usually evenly spread, are taken around the cargo or other objects being examined. However, this typically requires rotating the whole X-ray source and detector assembly as one integrated piece, around the cargo container, or rotating the cargo container inside a tunnel defined by X-ray source and detectors. Mechanical complexity, scanning speed, and other practical concerns impedes its adoption. In full scale CT, an object is viewed at multiple angles, at small increments, typically over at least 180 degrees. Analytic reconstruction techniques, such as filtered back projection, may be used to resolve information along X-ray beam paths, for example. Known full scale CT proposals have not been practical for scanning large objects because of their complexity.
Embodiments of the invention aim to provide improved imaging of objects, including large objects such as cargo containers. Other embodiments of the invention provide quasi-three-dimensional image reconstruction from viewer views than in CT reconstruction.
In one embodiment, accelerated electron beams are transported to multiple targets at respective locations to generate X-rays for radiography of an object at different viewing angles. In one example, resulting 3-D images enhance the ability to identify the structural contents of cargo containers. In addition, material signatures may be resolved along beam paths, greatly improving specificity by decreasing dilution due to overlapping, for example. Improved structural and material information (Z-value) improves automatic feature and threat detection.
In another embodiment, a system collimates radiation generated by a plurality of radiation sources into multiple beams for scanning an object at multiple angles. The object may be a cargo container, for example.
As discussed above, CT reconstruction algorithms, such as filtered backprojection (“FBP”), assume that data are continuously sampled along an ideal arc or line. In terms of physics, most also assume an ideal X-ray imaging mechanism, such as a noise-free monochromatic infinitesimal-width pencil-beam from a point source with no scatter using perfect noise-free linear electronics, for example. Both of these assumptions are problematic.
However, when the detector pixels are unlikely to all fall on a convenient curve or line, or to be finely spaced enough to be considered nearly continuous, or when the number of views is small or not finely spaced enough to be considered nearly continuous, analytic algorithms are generally precluded. Instead, iterative algorithms are typically used. However, in addition to geometry, it can also be very useful to include a more realistic physics model into the reconstruction process. While conventional CT often requires a host of correction algorithms (for beam hardening, cross-talk, noise reduction, scatter, non-linear gain, streak removal, etc), many of these corrections are ad hoc and do not completely solve the problem. Furthermore, many such artifacts come about largely due to the mismatch between the actual physics model and the ideal model assumed by the algorithm. Therefore, by using a more realistic physics model many of the artifacts may be reduced or eliminated, and the need for many of the usual corrections may be obviated.
Another class of algorithms, generally called Statistical Reconstruction (“SR”), allows for not only a more complete physics model, but also models the characteristics of the noise and uncertainty in the imaging process. SR also allows for some prior knowledge of the object, FBP and the original ART methods may be thought of as very particular statistical models in which the logarithm of the pixel values have are assumed to have uniform Gaussian noise, and all reconstructed objects are considered to be equally likely. However, neither assumption is particularly realistic, because X-rays contain polychromatic Poisson noise, and detectors add further Gaussian noise (on the detected signal, not its logarithm). This may be unimportant when a large number of X-rays pass through the object, but when few photons make it through the object, using an accurate noise model becomes important.
Examples of models for different materials and feasibility criteria to reduce candidate materials to yield useful material-identification information are described in U.S. Patent Publication No. 2013/0101156 A1, which is assigned to the assignee of the present invention and is incorporated by reference herein.
In accordance with an embodiment of the invention, an algorithm is described to generate a quasi-3D image from a set of projection measurements from a plurality of radiation sources in a plurality of different planes. Embodiments of the invention enable quasi-3D image reconstruction from much less information then required in CT. For example, conventional CT might use hundreds or thousands of views, whereas in accordance with embodiments of the invention, meaningful quasi-3D images may be obtained from the order of 5 to 20 views, for example. More views may be used, if desired.
In accordance with an embodiment of the invention, a method of reconstructing images is disclosed comprising scanning an object by a first plurality of radiation beams generated by a second plurality of sources at a plurality of angles and detecting the first plurality of radiation beams. Images of the object are reconstructed from projection measurements by searching for an estimate of an image that optimizes an imaging model defining a likelihood of the estimated image underlying the measured data, and a feasibility model defining a likelihood that the estimated image occurs in nature. The estimate of the image may comprise material information.
Multiple X-Ray Sources by Skewed Electron Transport
In the example of
The object 130 may be moved through the radiation scanning system 100 continuously or in steps, during scanning by a conveying system 150 indicated schematically by an arrow. The conveying system 150 may comprise a conveyor belt, a movable platform driven by rollers, or a towing system, for example, as is known in the art. Alternatively, the accelerator 110 and target stations 120a-120g may be mounted to a gantry and moved with respect to the object 130, continuously or in steps.
Three electron beam transport pipes 142, 144, 146 are provided in this example to convey the accelerated electrons to each of the target stations 120a-120g. A pipe 142 is shown at the near side of the object 130 and runs upward; another pipe 144, which is aligned with the path 112 of the electron beam from the accelerator 110, runs beneath object and across to the far side of the object; and a pipe 146 runs upward at the far side of the object. There are three target stations 120a, 120b, 120c along the pipe 142: station 120a near the top of the pipe 142 and object 130, station 120b near the middle of the pipe 142 and the object, and station 120c near the bottom of the pipe and the object.
There is one target station 120d along the pipe 144, near a middle of the pipe. There are three additional target stations 120e, 120f, 120g along pipe 146: target station 120e near the beginning of the pipe, target station 120f near the middle of the pipe, and target station 120g near a top and end of the pipe. The pipes 142-146 are skewed so that the positions of the seven X-ray target stations 120a-120g are offset, to scan the object 130 in different imaging planes. The radiation beams may be separated by one to two feet. This facilitates the placement of collimators (not shown) and detectors 190a-190g. Examples of detector configurations are shown in
An electromagnet 160 is provided between the pipes 142 and 144. A second electromagnet is 162 is provided between the pipes 144 and 146. Electromagnets 180a-180g are also provided at each of the seven X-ray target stations 120a-120f. Appropriate electromagnets are known in the art and are routinely used in medical accelerators, for example, to bend accelerated electrons. Such electromagnets are referred to as “bend magnets.” The electromagnets 160, 162, as well as the electromagnets 180a-180c, may be configured to bend the path of the electrons by 270°, for example. The electron beam pipes 142-146 are maintained at vacuum. Beam confinement and focusing can be applied along the pipes 142, 144, 146 to define a small electron beam spot as is also known in the art.
Detectors 190a-190g are provided on an opposite side of the object 130 as each target station 120a-120g, respectively, to detect radiation from each target station 120a-120g transmitted through the object 130, as discussed below and shown in
Magnetic fields are applied by the electromagnets 160, 162, to steer a pulse of accelerated electrons along a respective pipe 142, 146, toward the target stations along the respective pipes. When a proper magnetic field is applied by the electromagnet 160, the pulse of accelerated electrons is steered into the pipe 142. When no magnetic field is applied to the electromagnet 160, accelerated electrons go straight into the pipe 144 which is aligned with the electron beam path 112. To steer a pulse of accelerated electrons from the end of the pipe 144 (if coming this way) into the entrance end of the pipe 146, a proper magnetic field is applied by the electromagnet 162 to the accelerated electrons travelling through the pipe 144. A proper magnetic field is also always applied at the target station 120g to steer accelerated electrons to the target station (if coming this way).
When an interlaced dual energy accelerator 100 is used, as described in U.S. Pat. No. 8,183,801 B2, which is assigned to the assignee of the present invention and is incorporated by reference herein, for example, electron orbits through an orbit chamber may be kept the same for the high energy pulses and low energy pulses. The higher energy electrons need a higher magnetic field to bend at the same radius as lower energy electrons. Magnetic field strength is controlled by adjusting electric current through the electromagnet, as is known in the art. The electric current may be controlled by a processing device, such as the processing device 24 in
When a single energy accelerator 110 is used, accelerated electrons are cycled among the seven target stations 120a-120g. The following status Table A provides an example of magnet status for electron beam delivery to each target station 120a-120g. A blank means that its status does not matter, so the respective magnet can be on or off.
When an interlaced pulse dual-energy accelerator is used, the accelerator produces a pulse of higher energy electrons followed by a pulse of lower energy electrons, to alternately generate high energy radiation and lower energy radiation. For example, when scanning cargo containers, the high energy may be 9 MeV and the low energy 6 MeV or the high energy may be 6 MeV and the low energy may be 3.5 MeV. The required magnetic field is different for steering high energy electrons and low energy electrons. One possible pulsing sequence is:
As a result, a high energy radiograph and a low energy radiograph may be produced at each of the seven target stations 120a-120g. The number of radiographs doubles when two-layer stacked detectors 190a-190g are also used. Stacked detectors are described in U.S. Pat. No. 8,184,769, for example, which is assigned to the assignee of the present invention and is incorporated by reference herein. Radiographs from different X-ray target (source) stations 120a-120g have different viewing angles and may be used for image reconstruction, including tomographic reconstruction, to resolve information along X-ray beam paths (structural and material signature), as discussed further below.
While only one of the target stations 120a-120g is fired at a time, the pulsing sequence can run very quickly. For example when a 1,400 Hz accelerator is used, each X-ray target (source) station 120a-120g produces 100 high energy pulses and 100 low energy pulses per second.
As shown in
When an interlaced-pulse dual-energy accelerator 110 is used, the accelerator sequentially produces a pulse of higher energy electrons followed by a pulse of lower energy electrons at each target station 120a-120g. Slices of the object 130 are imaged at least twice at each X-ray target station 120a-120g, or from each viewing angle, at least once for each energy. All seven high energy radiographs are assembled into a high energy radiograph data set and all seven low energy radiographs are assembled into a low energy radiograph set. The number of radiograph sets doubles when two-layer stacked detectors are used, as discussed above. The interlaced accelerator may be a Mi6 Linatron® that generates radiation at 3.5 MeV and 6 MeV, or an Mi9 Linatron® that generates radiation at 6 MeV and 9 MeV, both available from Varian Medical Systems, Inc., Palo Alto Calif., for example. Three or more energies may be interlaced, as well. In addition, any number of layers may be provided in the stacked detectors.
Since the seven radiographs within any data set from different X-ray target stations 120a-120g have different viewing angles, they may be used for tomographic reconstruction to resolve structural and/or material information along the X-ray beam paths. One radiograph set (usually the high energy set for best quality) may be used to reconstruct a 3-D object distribution, slice by slice in this configuration, with a quasi 3-D reconstruction algorithm, for example. With structural information about the object 130, other radiograph sets may be used to further resolve material information (effective atomic number, Z value, or similar quantity) of each voxel. An example of a procedure quasi 3-D reconstruction algorithm for resolving this non-additive quantity is described below. Other reconstruction algorithms may be used, as well. While tomographic reconstruction may be provided, it is not required.
In-Plane Multiple X-Ray Sources by Electron Transport
Pipes 142, 144, 146 are provided, positioned adjacent to the same surfaces of the object 130, except that they lie in the same plane. The pipe 144 is aligned with the path 112 of the electrons output from the accelerator 110. The electromagnets 160 allows the accelerated electrons to proceed to the pipe 144 or directs the accelerated electrons to the pipe 142. The electromagnet 162 directs accelerated electrons through the pipe 146. Electromagnets 180a-180g direct accelerated electrons to the respective target stations 120a-120g. The detector array 302 comprises detector array sections 304-318, all lying in the same plane. Spaces are provided between adjacent detector array sections 304-318, respectively, to provide paths for the accelerated electrons bent by the electromagnets 180a-180g to the target stations 120a-120g, respectively, as shown in
At least a portion of all of the detector array sections 304-318 will detect radiation transmitted through the object 130 from each target station 120a-120g. Since only one target station 120a-120g works at one time, they are using the detectors in a time-sharing arrangement.
Integrating Out-of-Plane Scatter Detector
In one example, X-ray radiation generated from a respective target station 120a-120g is collimated into a fan beam. X-rays interact with the material in the object 130 along its beam path with a probability of interaction dependent on material type, material amount, and X-ray photon energy. Scattered X-ray photons travel along a path different from the original beam path. Most of scattered X-ray photons will travel out of the plane defined by the source collimator (not shown) and detectors 190a-190g, resulting in out-of-plane scatter. In addition, some X-ray photons will take a new path within the plane defined by source collimator and detectors 190a-190g resulting in in-plane scatter. While the mixture of scattered X-rays with non-attenuated X-rays (those that pass through the object 130 without interaction), may deteriorate imaging performance, scattered X-rays carry supplementary information about the object 130 that may be useful in image reconstruction if the scattered X-rays are recorded separately.
In accordance with an embodiment of the invention, out-of-plane scatter detector arrays provided adjacent to the transmission detector arrays 190a-190g to record out-of plane scattered X-rays, such as slightly out of plane scattered X-rays.
In one example, a proportion of the radiation recorded by the scatter detectors 402, 404, is subtracted from the radiation recorded by the respective detectors 190a-1190g, which detect radiation transmitted through the object 130.
Information from scatter detectors may be combined with that from transmission detectors and used in reconstruction algorithm by the processing device of the radiation scanning system 100, as discussed herein.
In-Plane Scatter Management with Rotating Detector Modulators
As mentioned above, with in-plane multiple X-ray sources, such as in the scanning system 300 of
Scatter detectors may be provided in this embodiment, in a similar manner as discussed above with respect to
The locations of target stations 120a-120g are examples and the stations may be located in other positions around the object 130 being scanned. In some examples, only two pipes may be provided.
While in these examples only one accelerator used to provide accelerated electrons to all the target stations 120a-120g, additional accelerators may be provided to provide accelerated electrons to at least two respective target stations. In this way, there are still fewer accelerators than target stations, saving costs.
Although straight detector arrays are shown in the description, curved detector arrays may be provided. Different pulse sequences (of energy and target station) other than the one described may also be used—as long as it covers all source stations and all energies in one cycle.
Quasi-3D Reconstruction
In accordance with another embodiment of the invention, quasi-3D image reconstruction is provided from a set of projection measurements resulting from radiation scanning of an object from a plurality of views by a respective plurality of radiation sources. These possibilities are referred to as quasi-3D imaging herein. The quasi-3-D reconstruction algorithms may be conducted by the image processor and the processing device of the scanning system, for example. The end result of a 3-D reconstruction algorithm shown to the user can be one of several forms:
Embodiments of the invention provide a more accurate geometric description of the scanning system, as well as a more realistic physics model, into the reconstruction process. Using feasibility information and prior knowledge of an object can also be useful. If the radiation scanning system itself does not generate much information, such as with an impenetrable X-ray path or a noisy scanner, there may be many images that are consistent with the scanner's projection data. However, many of those images might not be physically realistic. By narrowing choices down to those images that are physically likely, much less information is required from the scanning system. This observation has enabled new imaging methods, such as compressed sensing, and is also fundamental to recent advances in limited view CT reconstruction, especially for medical imaging in cases where human anatomy is well known and previous images of the same patient might also be available.
Model parameters related to the radiation scanning system are fit to the data, in Step 520. Model parameters may include the spectrum of the radiation beams and the geometry of the radiation scanning system, such as the geometrical relationship between the radiation sources, the detectors, and the object being examined, for example. The spectrum may be determined as described in U.S. Pat. No. 8,983,234, which was filed on Sep. 28, 2012, is assigned to the assignee of the present invention, and is incorporated by reference herein, for example. The geometrical relationship between the radiation source, the detectors, and the object 130 being examined may be determined as described in U.S. Pat. No. 7,922,390 B2 and U.S. Pat. No. 9,052,265, which are also assigned to the assignee of the present invention and are incorporated by reference herein, for example. The model parameters are stored in the calibration file, in Step 530.
The object is scanned, in Step 540, the calibration file is retrieved, in Step 550, and iterative reconstruction is performed in accordance with embodiments of the invention, in Step 560, as described below. A reconstructed image is provided to a viewer program, in Step 570, and scanning of the object continues in Step 540.
The optimization algorithm 730, in accordance with an embodiment of the invention, receives the initial guess. The optimization algorithm 730 searches for an underlying image that is most likely based on the imaging model 740 and the image feasibility model 750. The imaging model 740 describes how likely it is that measured data would have been collected, based on the current guess. The imaging model 740 is based, at least in part on a physics model related to all materials or to particular materials of interest, as discussed below. The imaging model 740 is also based on the calibration 760, which receives as inputs spectrum information 770 and the geometric information 780, described above. The image feasibility model 750 describes how likely it is that a given guess about an image can occur in nature, regardless of what is measured. These components are discussed in more detail below. The optimization algorithm 730 provides an output 790, which corresponds to Step 650 in
In security scanning, one might know many things about the objects being scanned. Some possibilities, referred to as feasibility criteria or “Regularization Options,” include:
In an example of an embodiment of the invention, a reconstruction algorithm performs the following steps:
1) Define P as the projection data acquired by a scanning system. Such as the scanning systems described herein or other scanning systems;
2) Define X as the quasi-3D image to be reconstructed;
3) Define D(X,P) as the data term, which is smaller when P and X are consistent, and larger when they are not;
4) Define R(X) as an unlikelihood term. It is smaller when X is reasonable, and larger when X is physically unlikely. When we are using several different types of prior knowledge, we may have several such functions R1, R2, . . . ;
5) Each unlikelihood term Rn is assigned some penalty weighting λn; and
6) Define C(X,P) as our reconstruction score, where C(X,P)=D(X,P)+ΣnλnRn(X). Then the reconstruction algorithm searches for the value of X that minimizes C(X,P).
Geometric Description
The geometric description discussed above with respect to Step 520 of
Typically, i=(x,y,z), where i is the coordinates of some reconstructed pixel or voxel X(i).
The amount of overlap of ray j and voxel I, may be described by v(i,j). Usually this is calculated as a volume, but it could alternately be a path length or a volume fraction.
Output Representation
The output of Step 650,
One possibility is to essentially make a separate quasi-3D reconstruction for each spectrum. That is Xk(i) is the product of (average volume density in voxel i) with (μkavg=average mass attenuation coefficient or the linear attenuation coefficient in spectrum k).
In this case, if the image is X then we would expect datum j to have a value similar to
where k is the effective spectrum used in ray j.
It may be useful to estimate density from X. Typically this is done by weighted average
where wk is some set of weights that sums to one.
This approach is also useful when there is only a single spectrum.
Material-Descriptors Representation
X may more directly characterize the materials. For megavolt cargo imaging the space of all possible materials appears to have a dimension around 3 to 5. If multiple spectra are used, K>5 may not help and K<3 may not capture all material information present in the imaging data. In one example K=4 may be used, for example.
A material model maps a material description to a corresponding interaction inside the physics model. In particular, it is useful for the material model to allow calculation of an approximate attenuation curve for a particular material. For a suitable set of basis functions, the approximate attenuation curve for voxel X is:
For a good set of basis functions and K≧4, in general for any real material, there is some choice of X(i) that makes the approximate attenuation curve very similar to the real attenuation curve. Such basis functions can be designed by a number of techniques. As described in U.S. Patent Publication No. 2013/0101156 A1, which is assigned to the assignee of the present invention and is incorporated by reference herein, one such technique is to form a matrix of mass attenuating values across all elements in the periodic table (or only elements of interest), in one matrix dimension and all energies of interest in another matrix dimension. A singular value decomposition (“SVB”) may then be performed on the matrices, and the largest singular vectors, such as the largest three singular vectors, the largest four singular vectors, or the largest five singular vectors, are kept. As is also described in U.S. Patent Publication No. 2013/0101156 A1, weighted SVD (“WSVD”) may also be performed, whereby energies are weighted by flux and detector sensitivity, materials are weighted by importance, and/or combinations of energy and material are weighted by expected transmission. WSVD is also described in “Weighted Low-Rank Approximations” by Srebo et al., in the Proceedings of the Twentieth International Conference on Machine Learning (ICML-2003), which is incorporated by reference herein.
A useful physics model (referred to in item 740 in
where φ(j,E) is the estimated spectrum for datum j. This can be known from the mechanical design, Monte Carlo simulation, or from a calibration method. Other Q functions can also be used (for example, if φ is too difficult to estimate, or the above model is not realistic enough). Other options include finding Q empirically as a function of X, during a calibration step, or using an analytic model (similar to the above) but modified to be either more realistic (perhaps more complicated) or to be easier to work with. In any case, the imaging model measures whether the measured data and the estimated data are both consistent with the physics model.
It is more complicated to estimate density from X in this version. Several options include:
ρ(i)=∥Xk(i)∥
ρ(i)=FXk(i) for some appropriately designed matrix F
The imaging model can be enacted via a data term that compares P and Q. Using likelihood estimation (regardless of the Q calculating method), the data term (Step 3, above), is:
which can be derived as the −log of the probability that we would have measured P given that we expected to measure Q and that our measurements are dominated by Poisson statistics. α is an optional weight that indicates the SNR of the spectrum associated with ray j (usually it is best to set this equal to the average effective number of photons through air, which can be estimated as α=SNR2 where SNR is measured in ray j through air).
There are also other alternatives for D. For example, a comparatively simple method is the simple least squares:
or least squares of logarithms:
For example, using least squares of logarithms, together with the “Separate-Spectrum” representation essentially leads to the basic ART algorithm.
Prior Knowledge
Regularization Option #1 can be expressed for separate-spectrum as
or for material-descriptors as
Regularization Option #2 can be expressed in a number of ways. Here is one possibility:
where ρ(i) is some estimate of density.
Regularization Option #3 can be expressed as
R(X)=Σ∥∇ρ(i)∥2 (easier, but gives blurry edges)
R(X)=Σ∥∇ρ(i)∥2 (more difficult, but gives sharper edges)
R(X)=Σψ(∇ρ(i)) (a compromise between the two)
Regularization Option #4 can be expressed by the analysis method using:
R(X)=∥Φtρ∥
where Φt is the analysis operator for an appropriate dictionary (such as wavelets).
One might also replace the norm operator with something more complicated such as a Huber function or entropy function.
Regularization Option #4 could instead be expressed by the synthesis method by, instead of representing X directly, storing a set of coefficients that can be used to synthesis X,
X=Φc
where Φ is some dictionary (such as wavelets or the like) that supports tensor fields.
Then we can use the penalty
R(X)=∥c∥
or we can replace the norm operator with any of the alternatives just mentioned.
In image processing, there are many known dictionaries Φ (whether used for analysis or synthesis), such as wavelets, curvelets, ridgelets, and the Fourier and discrete sine/cosine bases. There are also many known methods for designing dictionaries specifically tailored to particular types of images. Any such dictionaries can be used in embodiments of this invention (including building a dictionary specifically intended for cargo imaging).
Regularization Option #5 is an example of the burgeoning field of non-local image regularization, which is becoming popular for image restoration (denoising, or recovering missing pixels).
Regularization Option #6 can be expressed by
R(X)=Σ∥∇X∥2
R(X)=Σ∥∇X∥
or by summing the angles between neighbors, as in U.S. Patent Publication No. 2013/001156 A1, which is assigned to the assignee of the present invention, and is incorporated by reference herein. Any of the methods described in Holt, “Angular Regularization of Vector-Valued Signals”, ICASSP 2011, which is incorporated by reference herein, for example, may also be used.
Regularization Option #9 can be implemented at least in part by designing a matrix B then using the penalty
where details of how to construct such a B are also described U.S. Patent Publication No. 2013/0101156 A1, which is assigned to the assignee of the present invention and is incorporated by reference herein.
Optimization Algorithm
Once we have defined the D function and zero or more R functions, we can recover X by minimizing the cost C(X,P)=D(C,P)+ΣnλnRn as indicated in Step 6, above. A search algorithm known in the art may be used, or example. These include first-order gradient methods such as steepest descent, conjugate gradient search, or Nesterov's accelerated methods. It also could include second-order or quasi-second-order methods, including Newton-type methods and quasi-Newton methods (including L-BFGS). These also include proximal splitting methods, including augmented Lagrangian methods, ADMM and its variants, the Parallel-ProXimal Algorithm, the prox-lin method, iterative shrinkage, FISTA, NESTA, and the like. The choice of algorithm often depends on exactly what Regularizers are used (for example, if we switch between Regularization Option #3 vs #4 (analysis version) vs #4 (synthesis version), the fastest search algorithm might depend on which Regularization Option we choose).
There are also many prior art techniques for speeding up the search algorithm; these are generally complementary to the above techniques. These include using a direct/analytic algorithm as an initialization, multi-grid methods (i.e. coarse-to-fine progressive reconstruction), continuation methods (i.e. solve first for very large λ and gradually decrease), as well as hardware acceleration (using clusters, GPUs, FPGAs, ASICs, or the like).
Visualization & Output
If a 3D image is reconstructed (isotropic or anisotropic) for X, any of the volume rendering methods currently in use in medical or industrial imaging (including Multi-Planar Reformatting (MPR), Maximum-Intensity-Projection (MIP), or 3D ray-tracing-based algorithms (allowing a volume to either be first sliced, or to have a fractional transparency to each voxel), may be used.
If X is reconstructed as a 2D slice or a set of 2D slices, they may be viewed in a similar manner to viewing typical 2D images, except with additional user control over the slice location(s) and orientations. As the user changes the slice locations, the reconstruction can update the image. For 3D reconstruction, this full volume may be reconstructed and then the result may be manipulated (rotate, slice, etc) without re-reconstructing.
A radiography view and a material discrimination (“MD”) view may be provided. In the radiography view, a ρ image is calculated (using one of the X→ρ conversions listed above) and displayed to the user in grayscale. In the MD view, the X tensor is mapped into a more useful material space (such as atomic number) using any method from U.S. Patent Publication No. 2013/0101156 A1, which is assigned to the assignee of the present invention and is incorporated by reference herein. If there is confidence in the material assessment, the voxel may be displayed in an appropriate color, and if not, the voxel may be displayed in grayscale.
In addition to visualization, automated analysis may be performed on the images. Typically, different materials will have different X values. For example, if K=4, there is some 4-dimensional space that can be segmented into threat and non-threat regions (each of which might be made of many smaller disconnected regions). It can also be segmented directly into many smaller regions corresponding to more specific materials. In general these regions may overlap and may have soft boundaries. Voxels may be measured and classified as members of material-regions and colored appropriately (as a visualization technique) or a non-image-based response may be provided, such as an alarm bell, flashing icon, or text report, for example.
Another System for 2-D or Quasi-3D Cargo Imaging with Laminography or Dark-Field Imaging
In accordance with another embodiment of the invention, a radiation scanning system 800 (“system 800”) is disclosed comprising multiple sources and multiple detectors, to obtain multiple views in different planes, as shown in
On an opposite side of the object 814 are three detector arrays 818, 820, 822, with respective collimators 824, 826, 828, respectively.
The truck 816, and object 814 may be moved through the scanning system 800 by a conveyor system 830, as discussed above. Alternatively, radiation scanning system 800 may be supported by a gantry that moves across the object 814 and truck 816, as is known in the art.
The source collimators 808, 810, 812 in this example each define three (3) radiation beams at different angles. The detector arrays 818, 820, 822 are positioned in this example to record radiation transmitted through the object 814 from each radiation source 802, 804, 806. The collimators 808, 810, 812 define three (3) respective slots that are aligned with the radiation beams to be received by each detector array 818, 820, 822. Nine distinct views of the object 814 are thereby provided.
The radiation sources 802, 804, 806 and detector arrays 818, 820, 822 are not equally spaced so that all 9 rays are at different angles. The spacings between the radiation sources 802, 804, 806 and between the detector arrays 818, 820, 822 may be chosen to provide a fairly uniformly-spaced set of angles. The numbers of radiation sources 802, 804, 806 and the detector arrays 818, 820, 822 are just an example and different numbers of sources and detectors may be provided. For example, there may be one radiation source (emitting a wide beam) and at least two detector arrays to detect the wide beam, at least two radiation sources and one detector array, at least two radiation sources and at least three detector arrays, three radiation sources and four or five detector arrays, ten radiation sources and four detector arrays, twenty radiation sources and one detector array, or one radiation source and 20 detector arrays.
One of the source collimators 808, 810, 812 may define a radiation beam that is perpendicular to the truck's axis, to provide a direct side view. The radiation sources 802, 804, 806 may be all at the same elevation, so that the difference between images from each source/detector pair is only a rotation in the horizontal plane. Some of the radiation sources 802, 804, 806 may also be at different elevations, so that some images are effectively rotated axially with respect to each other.
Operation of the radiation sources 802, 804, 806 is timed by the processing device, such as the processing device shown in
Operation of the detectors 818, 820, 822 is coordinated with operation of the sources 802, 804, 806 by the processing device so that each time a source fires, a full measurement from each detector is collected. In the above example, each detector collects a full image every 1 msec.
Quasi-3D images may be reconstructed by the optimization algorithm, described above. 3D images may also be reconstructed from using computed laminography or tomosynthesis, or more generally, limited-view CT methods. Each view can also include multiple energies or stacked detectors, so we can incorporate material discrimination.
Scatter detectors may also be provided for dark-field imaging. Back-scatter detectors 830 and/or forward-scatter detectors 840 may be provided for example, as shown in
Measurements from the scatter detectors 830, 840, particularly the scatter detector 830, are more sensitive at the faces of the object 814 so they provide additional spatial information that can be incorporated into the same 3D reconstruction. A new image may be collected from each of the scatter detectors 830, 840 for each source pulse.
This embodiment of the invention uses cooperative imaging components, which may scale better than dual-view systems that use two non-interfering imaging chains, with one source and one detector each. Non-interfering designs scale linearly with the number of sources and detectors. In other words, in a conventional non-interfering design, buying N times as much imaging hardware gets you N times more data. However, if there are N sources in the present embodiment of the invention, adding only one more detector gives N times more information, and if there are M detectors, then adding only one more source gives M times more information. In other words, if M and N are kept similar, buying N times as much imaging hardware can get you N2 times more data, making it a far more cost effective option. In the example of
If, for example, to make useful volumetric reconstructions, a particular imaging task takes 25× more information than currently provided by a fourth generation scanner having a source rotating around an object, within a circular detector, then to scan a cargo container or other such object 814, in the same amount of time, a volumetric scanner using independent chains would be ˜25× more expensive than the current fourth generation scanner, which could be prohibitive. However, in embodiments of the present invention, the scanner would be ˜5× more expensive than the current fourth generation scanner, which is much more likely to become practically useful. Additional costs may be saved by sharing a power source, such as a modulator, and/or other components among the radiation sources 802, 804, 806.
Although example embodiments have been shown and described in this Specification and Figures, it would be appreciated by those skilled in the art that changes may be made to the illustrated and/or described example embodiments without departing from the scope of the invention, which is defined by the following claims.
The present application claims the benefit of U.S. Provisional Application No. 61/798,372, which was filed on Mar. 15, 2013, is assigned to the assignee of the present invention and is incorporated by reference herein. The present application is related to U.S. patent application Ser. No. 13/277,833, which was filed on Oct. 20, 2011 and was published on Apr. 25, 2013 bearing U.S. Patent Publication No. 2013/0101156, is assigned to the assignee of the present application is incorporated by reference herein.
Number | Name | Date | Kind |
---|---|---|---|
5504791 | Hell | Apr 1996 | A |
5654995 | Flohr | Aug 1997 | A |
5745546 | Hell | Apr 1998 | A |
5841046 | Rhodes et al. | Nov 1998 | A |
5841831 | Hell | Nov 1998 | A |
5940468 | Huang | Aug 1999 | A |
6035015 | Ruth | Mar 2000 | A |
6151381 | Grodzins | Nov 2000 | A |
6160869 | Zapalac | Dec 2000 | A |
6327339 | Chung | Dec 2001 | B1 |
6366641 | Whitham | Apr 2002 | B1 |
6369585 | Yao | Apr 2002 | B2 |
6407505 | Bertsche | Jun 2002 | B1 |
6421420 | Grodzins | Jul 2002 | B1 |
6493424 | Whitham | Dec 2002 | B2 |
6628745 | Annis | Sep 2003 | B1 |
6850592 | Schramm | Feb 2005 | B2 |
6937692 | Johnson | Aug 2005 | B2 |
7020241 | Beneke | Mar 2006 | B2 |
7060981 | Retterath | Jun 2006 | B2 |
7082182 | Zhou | Jul 2006 | B2 |
7099433 | Sommer | Aug 2006 | B2 |
7103137 | Seppi | Sep 2006 | B2 |
7203269 | Huber | Apr 2007 | B2 |
7208889 | Zavadtsev | Apr 2007 | B2 |
7208890 | Zavadtsev | Apr 2007 | B2 |
7215737 | Li | May 2007 | B2 |
7233644 | Bendahan | Jun 2007 | B1 |
7280631 | De Man | Oct 2007 | B2 |
7295651 | Delgado | Nov 2007 | B2 |
7319737 | Singh | Jan 2008 | B2 |
7333587 | De Man | Feb 2008 | B2 |
7356116 | Anwar | Apr 2008 | B2 |
7356118 | Might | Apr 2008 | B2 |
7379530 | Hoff | May 2008 | B2 |
7388940 | De Man | Jun 2008 | B1 |
7400701 | Cason | Jul 2008 | B1 |
7428292 | De Man | Sep 2008 | B2 |
7453987 | Richardson | Nov 2008 | B1 |
7505562 | Dinca | Mar 2009 | B2 |
7526064 | Akery | Apr 2009 | B2 |
7564943 | Sommer, Jr. | Jul 2009 | B2 |
7577236 | Nose | Aug 2009 | B2 |
7580505 | Kang | Aug 2009 | B2 |
7606348 | Foland | Oct 2009 | B2 |
7606349 | Oreper | Oct 2009 | B2 |
7616731 | Pack | Nov 2009 | B2 |
7627081 | Bontus | Dec 2009 | B2 |
7634045 | Popescu | Dec 2009 | B2 |
7639774 | De Man | Dec 2009 | B2 |
7646851 | Liu | Jan 2010 | B2 |
7684538 | Morton | Mar 2010 | B2 |
7688937 | Schomberg | Mar 2010 | B2 |
7702075 | Wang | Apr 2010 | B2 |
7706499 | Pack | Apr 2010 | B2 |
7724876 | Nose | May 2010 | B2 |
7813478 | Nisius | Oct 2010 | B2 |
7831012 | Foland | Nov 2010 | B2 |
7835486 | Basu | Nov 2010 | B2 |
7835499 | Yu | Nov 2010 | B2 |
7839967 | Grass | Nov 2010 | B2 |
7864917 | Ribbing | Jan 2011 | B2 |
7872241 | Rand | Jan 2011 | B2 |
7922390 | Holt et al. | Apr 2011 | B2 |
7991117 | Chen | Aug 2011 | B2 |
8000448 | Nose | Aug 2011 | B2 |
8137976 | Bjorkholm | Mar 2012 | B2 |
8138678 | Ishida | Mar 2012 | B2 |
8183801 | Chen | May 2012 | B2 |
8184769 | Fox | May 2012 | B2 |
8198587 | Whittum | Jun 2012 | B2 |
8203289 | Ho | Jun 2012 | B2 |
8222598 | Statham | Jul 2012 | B2 |
8223918 | Johnson | Jul 2012 | B2 |
8232748 | Treas | Jul 2012 | B2 |
8233682 | Fessler | Jul 2012 | B2 |
8284898 | Ho | Oct 2012 | B2 |
8300763 | Shedlock | Oct 2012 | B2 |
8311187 | Treas | Nov 2012 | B2 |
8422826 | Holt | Apr 2013 | B2 |
8457274 | Arodzero | Jun 2013 | B2 |
8508545 | Holt | Aug 2013 | B2 |
8509380 | Pelc | Aug 2013 | B2 |
8551785 | Bjorkholm | Oct 2013 | B2 |
8582720 | Morton | Nov 2013 | B2 |
8588494 | De Man | Nov 2013 | B2 |
8644549 | Foland | Feb 2014 | B2 |
8760050 | Mohr | Jun 2014 | B2 |
8761335 | Ishkhanov | Jun 2014 | B2 |
8774351 | Funk | Jul 2014 | B2 |
8803453 | Eaton | Aug 2014 | B2 |
8885975 | Yu | Nov 2014 | B2 |
8971484 | Beckmann | Mar 2015 | B2 |
9025024 | Xu | May 2015 | B2 |
9052265 | Holt | Jun 2015 | B2 |
9057679 | Morton | Jun 2015 | B2 |
9167681 | Cheung | Oct 2015 | B2 |
9257253 | Allen | Feb 2016 | B1 |
9258876 | Cheung | Feb 2016 | B2 |
9326366 | Krasnykh | Apr 2016 | B2 |
9330494 | Schultz | May 2016 | B2 |
20030194053 | Schramm et al. | Oct 2003 | A1 |
20070183568 | Kang et al. | Aug 2007 | A1 |
20070274456 | Holt | Nov 2007 | A1 |
20090041185 | Might et al. | Feb 2009 | A1 |
20100310042 | Fox et al. | Dec 2010 | A1 |
20100310175 | Holt | Dec 2010 | A1 |
20110096083 | Schultz | Apr 2011 | A1 |
20110116596 | Sommer, Jr. et al. | May 2011 | A1 |
20110129066 | Statham et al. | Jun 2011 | A1 |
20110280440 | Holt | Nov 2011 | A1 |
20120106816 | De Man et al. | May 2012 | A1 |
20130076913 | Xu et al. | Mar 2013 | A1 |
20130101156 | Holt | Apr 2013 | A1 |
Number | Date | Country |
---|---|---|
2004-174253 | Jun 2004 | JP |
2004-174261 | Jun 2004 | JP |
2009-261519 | Nov 2009 | JP |
WO-2008075595 | Jun 2008 | WO |
Entry |
---|
Fessler, Jeffrey A., “Iterative image reconstruction for CT”, EECS Dept., BME Dept., Dept. of Radiology University of Michigan, http://www.eecs.umich.edu/˜fessler, AAPM Image Educational Course—Image Reconstruction 11, Aug. 2, 2011. |
Holt, Kevin M., “Angular Regularization of Vector-Valued Signals”, ICASSP 2011. |
Srebro, Nathan et al., “Weighted Low-Rank Approximations”, Proceedings of the Twentieth International Conference on Machine Learning (ICML-2003). |
NPL Website, Published Aug. 31, 2007; http://www.classzone.com/books/earth—science/terc/content/investigations/esu101/esu101page05.cfm?chapter—no=investigation. |
Alvarez, Robert E. et al., “Energy-Selective Reconstructions in X-ray Computerized Tomography,” Phys. Med. Biol., vol. 21, No. 5, pp. 733-744; 1976. |
Office Action dated Apr. 19, 2017 issued in U.S. Appl. No. 13/277,833. |
Office Action dated Oct. 28, 2016 issued in U.S. Appl. No. 13/277,833. |
Office Action dated Mar. 4, 2015 issued in U.S. Appl. No. 13/277,833. |
Office Action dated Oct. 6, 2014 issued in U.S. Appl. No. 13/277,833. |
Advisory Action dated Jun. 23, 2014 issued in U.S. Appl. No. 13/277,833. |
Office Action dated Feb. 12, 2014 issued in U.S. Appl. No. 13/277,833. |
Office Action dated Sep. 18, 2013 issued in U.S. Appl. No. 13/277,833. |
Number | Date | Country | |
---|---|---|---|
20140294147 A1 | Oct 2014 | US |
Number | Date | Country | |
---|---|---|---|
61798372 | Mar 2013 | US |