Tomographic Imaging Using Fast Stochastic Hierarchical Bayesian Systems and Methods

0Associated
Cases 
0Associated
Defendants 
0Accused
Products 
0Forward
Citations 
0
Petitions 
1
Assignment
First Claim
1. An imaging system, comprising:
 a sensor configured to sense a scene from a plurality of different angles, thereby generating sensed data; and
an image generator device in communication with the sensor, wherein the image generator device is configured to;
generate a plurality of sensing matrices based on geometry of the angles and a waveform used to sense the scene,estimate a plurality of coefficients for an image based on the sensing matrices,estimate a Compound Gaussian (CG) vector based on the plurality of coefficients using a Simultaneous Perturbation Stochastic Approximation (SPSA), andgenerate an image based on the estimated CG vector.
1 Assignment
0 Petitions
Accused Products
Abstract
Systems and methods are provided for imaging that demonstrably outperform previous approaches (especially compressive sensing based approaches). Embodiments of the present disclosure provide and solve an imaging cost function via a stochastic approximation approach. By doing so, embodiments of the preset disclosure provide a significant means of generalization and flexibility to adapt to different application domains while being competitive in terms of computational complexity.
0 Citations
No References
No References
20 Claims
 1. An imaging system, comprising:
a sensor configured to sense a scene from a plurality of different angles, thereby generating sensed data; and an image generator device in communication with the sensor, wherein the image generator device is configured to; generate a plurality of sensing matrices based on geometry of the angles and a waveform used to sense the scene, estimate a plurality of coefficients for an image based on the sensing matrices, estimate a Compound Gaussian (CG) vector based on the plurality of coefficients using a Simultaneous Perturbation Stochastic Approximation (SPSA), and generate an image based on the estimated CG vector.  View Dependent Claims (2, 3, 4, 5, 6, 7)
 8. A method for generating an image, the method comprising:
generating, using an image generator device, a plurality of sensing matrices based on geometry of a plurality of angles used to sense a scene and a waveform used to sense the scene; estimating, using the image generator device and based on the sensing matrices, a plurality of coefficients for an image; estimating, using the image generator device, a Compound Gaussian (CG) vector based on the plurality of coefficients using a Simultaneous Perturbation Stochastic Approximation (SPSA); and generating, using the image generator device, the image based on the estimated CG vector.  View Dependent Claims (9, 10, 11, 12, 13, 14, 15)
 16. An imaging device, comprising:
a processor; a memory; and an image generator device configured to generate a plurality of sensing matrices based on geometry of a plurality of angles used to sense a scene and a waveform used to sense the scene, estimate, based on the sensing matrices, a plurality of coefficients for an image, estimate a Compound Gaussian (CG) vector based on the plurality of coefficients using a Simultaneous Perturbation Stochastic Approximation (SPSA), and generate an image based on the estimated CG vector.  View Dependent Claims (17, 18, 19, 20)
1 Specification
This application claims the benefit of U.S. Provisional Patent Application No. 62/751,742, filed on Oct. 29, 2018, which is incorporated by reference herein in its entirety.
The United States Government has ownership rights in this invention. Licensing inquiries may be directed to Office of Technology Transfer at US Naval Research Laboratory, Code 1004, Washington, D.C. 20375, USA; +1.202.767.7230; techtran@nrl.navy.mil, referencing Navy Case Number 107290US3.
This disclosure relates to imaging, including tomographic imaging.
A wide variety of analytical and statistical approaches have proliferated over the past several decades in diverse imaging applications, such as radar, medical, and astronomical imaging. Image recovery algorithms attempt to achieve the highest quality reconstruction in a timely manner. Robust and highfidelity image reconstruction can be difficult, especially given a limited number of samples and in the presence of sensor noise. What is needed are systems and methods for imaging that take dramatically fewer operations while retaining high reconstruction quality.
The accompanying drawings, which are incorporated in and constitute part of the specification, illustrate embodiments of the disclosure and, together with the general description given above and the detailed descriptions of embodiments given below, serve to explain the principles of the present disclosure. In the drawings:
Stochastic Hierarchical Bayesian MAP (fsHBMAP) algorithm in accordance with an embodiment of the present disclosure;
Features and advantages of the present disclosure will become more apparent from the detailed description set forth below when taken in conjunction with the drawings, in which like reference characters identify corresponding elements throughout. In the drawings, like reference numbers generally indicate identical, functionally similar, and/or structurally similar elements. The drawing in which an element first appears is indicated by the leftmost digit(s) in the corresponding reference number.
In the following description, numerous specific details are set forth to provide a thorough understanding of the disclosure. However, it will be apparent to those skilled in the art that the disclosure, including structures, systems, and methods, may be practiced without these specific details. The description and representation herein are the common means used by those experienced or skilled in the art to most effectively convey the substance of their work to others skilled in the art. In other instances, wellknown methods, procedures, components, and circuitry have not been described in detail to avoid unnecessarily obscuring aspects of the disclosure.
References in the specification to “one embodiment,” “an embodiment,” “an exemplary embodiment,” etc., indicate that the embodiment described may include a particular feature, structure, or characteristic, but every embodiment may not necessarily include the particular feature, structure, or characteristic. Moreover, such phrases are not necessarily referring to the same embodiment. Further, when a particular feature, structure, or characteristic is described in connection with an embodiment, it is submitted that it is within the knowledge of one skilled in the art to understand that such description(s) can affect such feature, structure, or characteristic in connection with other embodiments whether or not explicitly described.
Embodiments of the present disclosure provide systems and methods for imaging that demonstrably outperform previous approaches (especially compressive sensing based approaches). Embodiments of the present disclosure provide and solve an imaging cost function via a stochastic approximation approach. By doing so, embodiments of the present disclosure offer a significant means of generalization and flexibility to adapt to different application domains, while being competitive in terms of computational complexity.
For example, embodiments of the present disclosure systematically incorporate prior (probabilistic) knowledge of the scene structure into the inference process. Embodiments of the present disclosure use a robust statistical inference technique under subNyquist sampling when there is degradation of the received signal due to corruptions arising from environmental and other factors.
Image recovery algorithms attempt to achieve the highest quality reconstruction in a timely manner. The former can be achieved in several ways, among which are by incorporating Bayesian priors that exploit natural image tendencies to cue in on relevant phenomena. The Hierarchical Bayesian Maximum a Posteriori (HBMAP) is one such approach which produces compelling results albeit at a substantial computational cost. Embodiments of the present disclosure largely retain the imaging quality of HBMAP while achieving an order of magnitude reduction in computational speed via the application of stochastic approximation techniques.
For example, embodiments of the present disclosure retain the proficient nature of Hierarchical Bayesian Maximum a Posteriori'"'"'s (HBMAP) TypeI estimation and use a stochastic approximationbased approach to TypeII estimation. The resulting algorithm, Fast Stochastic Hierarchical Bayesian MAP (fsHBMAP), takes dramatically fewer operations while retaining high reconstruction quality. Embodiments of the present disclosure employ a fsHBMAP scheme to the problem of tomographic imaging and demonstrate that fsHBMAP furnishes promising results when compared to many competing methods.
In an embodiment, an fsHBMAP algorithm is used to provide a vehicle to extend the capabilities of imaging algorithms to include truly global priors. Probabilistic graphical models can yield stateoftheart performance for signal classification problems. Embodiments of the present disclosure deploy CG based probabilistic graphical models for the purposes of regression to obtain stateoftheart results. The formulation and fsHBMAP algorithm used by embodiments of the present disclosure can be applied to a variety of sampling matrices and can thus be potentially applied to a wide range of inverse problems beyond imaging.
The act of obtaining images from tomographic settings has been a long studied problem with applications ranging from radar to magnetic resonance imaging. However the signaling apparatus is constructed, there is a consistent theme: a plane/volume is captured by rays of directed waves that are recorded, match filtered, and processed to achieve the pixels of the image. Several aspects to this procedure have emerged as especially difficult, none the least of which is limited sampling. For example, one problem is obtaining a quality image if a tomographic setup can collect only a dearth of jittered data points.
3. Exemplary fsHBMAP Approach
Generally, a received signal is equal to the generalized Radon transform of an image plus noise. An unknown image (I) can be represented, in terms of basis functions, as I=Φc, where Φ∈^{p×n }is the Dictionary matrix in which the image is represented (e.g., for Wavelets etc.). In an embodiment, Φ can be chosen to be the Wavelet dictionary, and, given this, the only unknown quantities to be estimated are the coefficients c∈^{n}. However, it should be understood that 0 is not limited to wavelets in accordance with embodiments of the present disclosure.
Given the inputs discussed above with reference to
In an embodiment, based on the above, the wavelet coefficients of an image (e.g., subimage 206) can be calculated.
In an embodiment, the Compound Gaussian (CG) Model of c can be represented by c=z⊙u, where u˜N(0, σ_{u}) is a Gaussian random vector, and z≥0 is a NonGaussian random vector. In an embodiment, vector z can be estimated (given c in this case), and vector c can be normalized (u(i)=c(i)/z(i) for each vector component i). Next, in an embodiment, the histogram of u and its best Gaussian match can be plotted.
In an embodiment, Simultaneous Perturbation Stochastic Approximation (SPSA) approach can be used for stochastic approximation. In an embodiment, SPSA (and related approaches such as Finite Difference SA (FDSA)) starts with finite difference approximations (e.g., of the gradient) and incorporates probabilistic theories to reduce computations. Embodiments of the present disclosure can drastically speed up the estimation of the nonGaussian component (i.e., the zcomponent) via the application of SPSA.
While fsHBMAP techniques are discussed above with reference to image reconstruction, it should be understood that fsHBMAP techniques can be used to generate original image(s) based on information from sensors. For example, in an embodiment, a scene/environment can be sensed from a variety of different angles (e.g., using sensors on board a ship, such as radar, syntheticaperture radar (SAR), etc.), and fsHBMAP techniques can be used to generate an original image based on information from these sensors.
In an embodiment, image generator 706 (or, in an embodiment, sensor(s) 702) generates a plurality of sensing matrices (e.g., in an embodiment, using a Radon transform) based on based on geometry of the sensing angles and the sensor waveform used to sense the scene. In an embodiment, image generator 706 estimates a plurality of coefficients for an image based on the sensing matrices. In an embodiment, image generator 706 estimates a Compound Gaussian (CG) vector based on the plurality of coefficients using a Simultaneous Perturbation Stochastic Approximation (SPSA) and generates an image 712 based on the estimated CG vector.
In an embodiment, a Compound Gaussian (CG) model c of the coefficients can be represented by the equation c=z⊙u, where z is the nonGaussian component and u is the Gaussian component. Thus, in an embodiment, both the Gaussian component u and the nonGaussian component z are estimated to generate an estimated CG model c of the coefficients.
For example, in an embodiment, image generator 706 estimates a nonGaussian component z of a Compound Gaussian (CG) vector c based on the plurality of coefficients using a Simultaneous Perturbation Stochastic Approximation (SPSA). In an embodiment, image generator 706 multiplies columns of an effective sensing matrix with scalar components of the nonGaussian component z of the CG vector c. For example, in an embodiment, an effective sensing matrix is generated from the matrix product of the sensing matrix and a Dictionary matrix. As discussed above, the sensing matrix can be generated based on a Radon transform, and Φ can be chosen to be the Wavelet dictionary.
In an embodiment, image generator 706 estimates a Gaussian component u of the CG vector c based on the normalized effective sensing matrix (e.g., after it was multipled), and the sensed data. In an embodiment, image generator 706 estimates a best Gaussian component u of the CG vector c. For example, in an embodiment, the best Gaussian component u is selected as the least square estimate of the Gaussian component of the CG vector c. In an embodiment, image generator 706 can then generate the image 712 based on the estimated CG vector c and a dictionary matrix.
Imaging device 704 and/or image generator 706 can be implemented using hardware, software, and/or a combination of hardware or software. Imaging device 704 and/or image generator 706 can be implemented using a single device or using multiple devices. Imaging device 704 and/or image generator 706 can be implemented as standalone (specialpurpose) devices or can be integrated into a host device. In an embodiment, sensor(s) 702 can be integrated into imaging device 704. In an embodiment, sensor(s) 702 can be integrated into image generator 706. In an embodiment, image generator 706 is configured to receive data from a source other than sensor(s) 702. For example, in an embodiment, image generator 706 can be configured to receive data for image processing from a computer, from a network (e.g., the Internet), from a wireless link (e.g., Bluetooth), or from any other source.
In step 808, a nonGaussian component of a Compound Gaussian (CG) vector is estimated (e.g., by image generator 706) based on the plurality of coefficients using a Simultaneous Perturbation Stochastic Approximation (SPSA). In step 810, columns of an effective sensing matrix are multiplied (e.g., by image generator 706) with scalar components of the nonGaussian component of the CG vector. In step 812 a Gaussian component of the CG vector is estimated (e.g., by image generator 706) based on the effective sensing matrix, the nonGaussian component, and the sensed data. In step 814, an image is generated (e.g., by image generator 706) based on the estimated CG vector and a dictionary matrix.
Exemplary algorithms for HBMAP and fsHBMAP will now be discussed in greater detail. Of the different manners in which to handle such a situation, the Hierarchical Bayesian Maximum a Posteriori (HBMAP) approach has yielded compelling results on optical images with the express intent of usage for, but not limited to, tomography. However this method has a significant drawback—it is computationally expensive. Its detailed structure involves a complicated objective function which requires nontrivial computational cost, restricting its applicability to largescale inverse problems. A brief overview of an exemplary HBMAP algorithm will now be discussed.
In an embodiment, the following linear model can be considered for reconstructing images, which finds powerful applications in imaging and compressive sensing: y=Xβ+ϵ, where y is a vector of received measurements, X is a dictionary containing some basis, whether it be wavelet, Fourier, or similar, β is the coefficient vector we look to find, and ϵ is Gaussian white noise. In an embodiment, we presume y, β, ϵ∈^{n }and X∈^{n×n}. There is no shortage in ways to solve for β given the assumption of the above linear model, and many involve a prior being placed on this coefficient term either explicitly or implicitly via the inclusion of a penalty function. In this way, it can be further interpreted as solving Equation (1) below:
In an embodiment, a good choice for the prior on β for this distribution is the Laplacian, i.e., logP(β)∝−β_{1}, which can be interpreted as a relaxation of sparsityenforcing to pseudo norm penalty. This choice for β lacks proper sophistication in order to capture certain natural phenomena and, in an embodiment, is where the HBMAP comes into play. First, let us revisit a linear model construction. In an embodiment, for HBMAP, X can be thought of as the product of two quantities: a measurement matrix ψ (e.g., a sensing matrix) and a dictionary Φ. In an embodiment, with X=ψΦ, β is a coefficient vector (with the true image being Φβ). This grants the flexibility to further incorporate convenient linear transformation spaces (typically, but not restricted to, wavelet) concerning the image I. In an embodiment, ψ could hold the information regarding the arrangement of sensors and readings, primarily in the mold of a Radon transform matrix in the tomographic setting, and Φ a wavelet dictionary such that Φβ is the optical image to be recovered. In an embodiment, a Compound Gaussian (CG) can be a suitable prior for β since it is now the wavelet coefficients of a natural image. In an embodiment, β can be shown by Equation (2) below:
β=u⊙z where u˜(0, Σ_{u}) and z=h(x) (2)
In Equation (2), x follows some mutliscale Gaussian tree structure, and h(⋅) is an entrywise, usually nonlinear function. In an embodiment, this construction leads to a two stage solution starting with a calculation for z (TypeII estimation), followed by a scheme for u (TypeI estimation). In an embodiment, the TypeII estimation involves substantial computations. In an embodiment, the issues arise from x; we can state y as shown by Equation (3) below:
y=X(h(x)⊙u)+ϵ (3)
In an embodiment, optimizing z is the same as optimizing x as the input to the nonlinear h. In an embodiment, TypeII estimation, due to the multiscale Gaussian hierarchical Bayesian structure imposed on x, comes down to solving Equations (4) and (5) below given H(x)=diag(h(x)):
In an embodiment, this objective function, which we will define as f(x), has what could be described as a particularly nasty gradient in terms of computational efficiency. Even with a “nice” choice for h (in terms analytical properties such as differentiability), the gradient can require costly Kronecker products, which are unavoidable without severe modification. Embodiments of the present disclosure can circumvent the calculation of an unseemly ∇_{x}f and an even more complex ∇_{x}^{2}f.
For example, assuming that an optimal x* is obtained, TypeI estimation can build off of this to calculate u. In an embodiment, this in tum is much simpler, as shown by Equation (6) below:
In an embodiment, for Λ=diag(z*), Equation (6) leads to Equation (7) below:
Λ(X^{T}Σ_{ϵ}^{−1}X+Λ^{−1}Σ_{u}Λ^{−1})Λu=ΛX^{T}Σ_{ϵ}^{−1}y (7)
In an embodiment, since this computation for u can be nontrivial in the manner given by (7), a strategy can be employed to reduce the load by defining Λ_{τ }as shown by Equation (8) below:
In an embodiment, based on the above, u can be found using Equation (9) below:
Λ_{96 }(X^{T}Σ_{ϵ}^{−1}X+Λ^{−1}Σ_{u}Λ^{−1})Λ_{τ}u=ΛX^{T}Σ_{ϵ}^{−1}y (9)
In an embodiment, with both z=h(x*) and u, one can quickly find β using Equation (2), completing the two step process.
7. Exemplary fsHBMAP Algorithm
Embodiments of the present disclosure provide a modified HBMAP algorithm—a fast stochastic HBMAP (fsHBMAP), using a stochastic approximation approach, for solving a problematic objective function. An exemplary fsHBMAP algorithm in accordance with an embodiment of the present disclosure will now be discussed with reference to the HPMAP algorithm discussed above.
In an embodiment, in Equation (7), if no simplifying substitutions regarding Λ are made, we the problem can be reduced to Equation (10) below:
In Equation (10), Σ_{ϵ}=σ_{ϵ}^{2}I, and Σ_{u}=σ_{u}^{2}I. Notice that, in Equations (11) and (12):
In an embodiment, the optimality of β hinges on the assumed optimality of x* (i.e., where h(x*) constitutes the diagonal of Λ. In an embodiment, a suboptimal x can be used as the input to h. In an embodiment, if the value of σ_{ϵ}^{2 }is known, or, more likely, estimated well, then {tilde over (β)} the coefficient vector corresponding to the solution of Equation (12) using a suboptimal x, can perform nontrivially.
In an embodiment, given this aspect of HBMAP, the question then becomes what influence the TypeII estimation has over the algorithm'"'"'s output. Referencing
As mentioned above, TypeII estimation involves a prohibitively costly gradient calculation (to say nothing of its Hessian). Embodiments of the present disclosure provide an approximation {tilde over (f)} that is much more efficient and still impressively capable. To start, let us break up f, as shown by Equation (13):
As can be seen by inspection, f_{3 }represents a more normallooking objective function. In an embodiment, key technical issues revolve around ameliorating the numerical tractability of f_{1 }and f_{2 }given their role in B(x) as defined in Equation (5). Regarding B(x), the nonlinear choice of h can be used, as shown by Equation (14):
In Equation (14), a∈ (note that a framework in accordance with embodiments of the present disclosure can use a general choice for both). Substituting this and our other assumptions for the matrices Σ_{u }and Σ_{ϵ }into Equation (5) yields Equation (15) below:
B(x)=σ_{u}^{2}X^{T}H^{2}(x)X+σ_{ϵ}^{2}I (15)
In an embodiment, a slightly simplified B(x) is introduced, and a formal choice for H is performed for f_{2}. In an embodiment, two more assumptions are made: ψ is square and full rank. The proceeding theory and experimental work provides a strong basis for these assumptions. By the Minkowski inequality, as shown by Equation (16):
In an embodiment, for the determinant of the matrix X, the dictionary D is typically chosen so that it is unitary, meaning that det(D)=1. Thus, it can be seen that det(X)=det(RD)=det(R). Therefore, since log is monotonically increasing, this leads to Equation (17) below:
det(σ_{u}^{2}X^{T}H^{2}(x)X)=σ_{u}^{2n}det(R)^{2}det(H^{2}(x)) (17)
Which, returning to the log det function, leads to Equation (18) below:
In Equation (18), K=log(σ_{u}^{2n}det(R)^{2})+log(σ_{ϵ}^{2n}). In an embodiment, an approximation for f_{2 }can be defined as
In an embodiment, moving onto f_{1 }reveals no readily actionable approximation. The inverse of B(x), which is composed of the sum of matrices, is a highly nontrivial problem. Embodiments of the present disclosure avoid gradient calculations altogether by invoking a stochastic approximation (SA) approach.
In an embodiment, if we were to start with the idea that we just wanted a near enough calculation of the gradient for use in a steepest descent technique, then it would not be unreasonable to pursue some finite difference scheme, and SA builds off this idea. In an embodiment, for SPSA, the update step of a gradient descent algorithm is modified for step size α_{k}, as shown by Equation (19) below:
x_{k+1}=x_{k}−α_{k}{tilde over (g)}(x) (19)
In Equation (19), there is a gradient approximation, {tilde over (g)}, as shown by Equation (20) below:
In Equation (20), c_{k}>0 is a small value that decreases as k→∞, Δ_{k }follows a symmetric±1 Bernoulli distribution, and is the elementwise division operator. Given the smoothness off and sufficient choices for α_{k }and c_{k}, we found that SPSA served as competent method for TypeII estimation with substantially fewer calculations. Note that SPSA involves only two evaluations of the loss function at each iteration.
To gain an idea of how helpful SPSA can be and what computation effort can be expected, refer to
To provide context for fsHBMAP, we experimented on a 32×32 piece 504 of image 202 projected with a sampling pattern containing 18 rays with 71 samples each. The result 504 obtained fsHBMAP had a SSIM of 0.78. The result 504 from fsHBMAP in accordance with an embodiment of the present disclosure outperforms competing methods with a substantially less computation time. Further, embodiments of the present disclosure using fsHBMAP are much more amenable to being applied to large scale imaging problems because, in an embodiment, no explicit gradient or Kronecker calculations are involved.
It is to be appreciated that the Detailed Description, and not the Abstract, is intended to be used to interpret the claims. The Abstract may set forth one or more but not all exemplary embodiments of the present disclosure as contemplated by the inventor(s), and thus, is not intended to limit the present disclosure and the appended claims in any way.
The present disclosure has been described above with the aid of functional building blocks illustrating the implementation of specified functions and relationships thereof. The boundaries of these functional building blocks have been arbitrarily defined herein for the convenience of the description. Alternate boundaries can be defined so long as the specified functions and relationships thereof are appropriately performed.
The foregoing description of the specific embodiments will so fully reveal the general nature of the disclosure that others can, by applying knowledge within the skill of the art, readily modify and/or adapt for various applications such specific embodiments, without undue experimentation, without departing from the general concept of the present disclosure. Therefore, such adaptations and modifications are intended to be within the meaning and range of equivalents of the disclosed embodiments, based on the teaching and guidance presented herein. It is to be understood that the phraseology or terminology herein is for the purpose of description and not of limitation, such that the terminology or phraseology of the present specification is to be interpreted by the skilled artisan in light of the teachings and guidance.
Any representative signal processing functions described herein can be implemented using computer processors, computer logic, application specific integrated circuits (ASIC), digital signal processors, etc., as will be understood by those skilled in the art based on the discussion given herein. Accordingly, any processor that performs the signal processing functions described herein is within the scope and spirit of the present disclosure.
The above systems and methods may be implemented using a computer program executing on a machine, a computer program product, or as a tangible and/or nontransitory computerreadable medium having stored instructions. For example, the functions described herein could be embodied by computer program instructions that are executed by a computer processor or any one of the hardware devices listed above. The computer program instructions cause the processor to perform the signal processing functions described herein. The computer program instructions (e.g., software) can be stored in a tangible nontransitory computer usable medium, computer program medium, or any storage medium that can be accessed by a computer or processor. Such media include a memory device such as a RAM or ROM, or other type of computer storage medium such as a computer disk or CD ROM. Accordingly, any tangible nontransitory computer storage medium having computer program code that cause a processor to perform the signal processing functions described herein are within the scope and spirit of the present disclosure.
While various embodiments of the present disclosure have been described above, it should be understood that they have been presented by way of example only, and not limitation. It will be apparent to persons skilled in the relevant art that various changes in form and detail can be made therein without departing from the spirit and scope of the disclosure. Thus, the breadth and scope of the present disclosure should not be limited by any of the abovedescribed exemplary embodiments.