Next Article in Journal
Combination of High- and Low-Rate GPS Receivers for Monitoring Wind-Induced Response of Tall Buildings
Next Article in Special Issue
Applicability of Diffuse Ultrasound to Evaluation of the Water Permeability and Chloride Ion Penetrability of Cracked Concrete
Previous Article in Journal
Location-Based Lattice Mobility Model for Wireless Sensor Networks
Previous Article in Special Issue
Improvement of Detection Sensitivity of Microbubbles as Sensors to Detect Ambient Pressure
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Sparse Ultrasound Imaging via Manifold Low-Rank Approximation and Non-Convex Greedy Pursuit

by
Thiago Alberto Rigo Passarin
*,†,
Marcelo Victor Wüst Zibetti
and
Daniel Rodrigues Pipa
Graduate Program in Electrical and Computer Engineering (CPGEI), Federal University of Technology, Paraná (UTFPR), Curitiba PR 80230-901, Brazil
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Sensors 2018, 18(12), 4097; https://doi.org/10.3390/s18124097
Submission received: 19 October 2018 / Revised: 16 November 2018 / Accepted: 20 November 2018 / Published: 23 November 2018
(This article belongs to the Special Issue Ultrasonic Sensors 2018)

Abstract

:
Model-based image reconstruction has improved contrast and spatial resolution in imaging applications such as magnetic resonance imaging and emission computed tomography. However, these methods have not succeeded in pulse-echo applications like ultrasound imaging due to the typical assumption of a finite grid of possible scatterer locations in a medium–an assumption that does not reflect the continuous nature of real world objects and creates a problem known as off-grid deviation. To cope with this problem, we present a method of dictionary expansion and constrained reconstruction that approximates the continuous manifold of all possible scatterer locations within a region of interest. The expanded dictionary is created using a highly coherent sampling of the region of interest, followed by a rank reduction procedure. We develop a greedy algorithm, based on the Orthogonal Matching Pursuit, that uses a correlation-based non-convex constraint set that allows for the division of the region of interest into cells of any size. To evaluate the performance of the method, we present results of two-dimensional ultrasound imaging with simulated data in a nondestructive testing application. Our method succeeds in the reconstructions of sparse images from noisy measurements, providing higher accuracy than previous approaches based on regular discrete models.

1. Introduction

Model-based image reconstruction methods provided important advances to imaging techniques such as magnetic resonance imaging (MRI) [1] and emission computed tomography (ECT) [2] in the last few decades. These methods rely on a known model that results in the captured signal being represented by a sum of N coefficient-weighted responses. These responses are usually point spread functions (PSF), and coefficients are usually intensities of pixels at the modelled locations. Then, by using the discrete model, a vector of acquired data and a regression algorithm, the intensity of each pixel is determined [3]. The use of model-based techniques in ultrasound imaging relies on the assumption that all reflectors (or scatterers) are located on any of a finite grid of N modelled positions [4].
Real-world inspected objects easily break this assumption and many scatterers may be situated off-grid. Many previous studies with model-based algorithms for ultrasound imaging [4,5,6,7,8,9,10,11] have reported that resolution and contrast are substantially improved in comparison to delay-and-sum (DAS) algorithms when data comes from simulations with scatterers located strictly on a modelled grid. However, images are corrupted by artifacts when the grid is not respected, which is typical in data acquired from real measurements. Consequently, DAS beamforming algorithms remain as state-of-the-art for ultrasound imaging, despite having well understood physical limitations regarding spatial resolution [12,13].

2. Model-Based Imaging and Regularization

Let R M be the space of the data observed through an acquisition process. A single, unity amplitude event located at position τ R D causes the discrete acquired signal y ( τ ) R M , known as the PSF. In ultrasound imaging, the event denotes a point-like reflexivity (also called a scatterer) [14,15], as represented in Figure 1, and D typically equals 2 as the reflexivity is being mapped over a two-dimensional plane. The variation of the set of D parameters τ within a region of interest describes a D-dimensional manifold
M : = { y ( τ ) : τ ROI }
of all possible PSFs on R M . We will consider the two-dimensional case, where τ = [ x , z ] T (being · T the transpose) are the lateral and axial spatial dimensions, respectively.
An acquired signal c R M is assumed to be composed by a sum of N events, or N samples, from the continuous PSF manifold
c = n = 1 N v n y ( x n , z n ) + w ,
where v n is the amplitude of the n-th event and the vector w R M accounts for acquisition noise, which we will assume to be Gaussian white noise with variance σ 2 .
In a pulse-echo image with N pixels, v n in Equation (2) encodes the reflexivity of the n-th scatterer, located at position ( x n , z n ) , and is represented as the brightness of the corresponding pixel. This implies sampling the parameters ( x , z ) on N possible scatterer locations (or pixels).
Considering N coordinate pairs ( x n , z n ) , we make h n = y ( x n , z n ) , n = 1 , , N , and define the model matrix H = [ h 1 , , h N ] R M × N . Then, Equation (2) can be written in compact form as
c = H v + w ,
where v = [ v 1 , , v N ] T is the vector of scatterer amplitudes. This model has been used in B-mode (two-dimensional) [4,5,6,7,8,9], A-mode (one-dimensional) [16,17], and three-dimensional [18] ultrasound imaging.
The reconstruction of the amplitudes vector v from a given acquisition c in Equation (3) is based on the minimization of a cost function, such as the least squares (LS) problem
v ^ = arg min v c H v 2 2 ,
which is linear and can be solved by well-known methods [19].
However, real-world matrices are often ill-conditioned, which causes artifacts on the reconstructed images in the presence of noise [20], even when all events are on grid. In ultrasound imaging, this problem has been addressed with linear regularization methods such as Truncated SVD (TSVD) [7] and Tikhonov regularizarion [5,6,8].
Sparsity-promoting regularization penalties such as p -(pseudo)norm minimization with p 1 have shown successful results in ultrasound non-destructive testing (NDT), where the assumption of sparsity in the space domain reflects the nature of discontinuities in observed materials [4,9,17,21].
Greedy algorithms effectively solve reconstruction problems where the cost function involves the 0 pseudonorm. In [10], sparsity is induced in the solution by the assumption that the presence of scatterers can be modelled by a Bernoulli process with a low value for the probability parameter. The problem is then solved with a greedy algorithm called Multiple Most Likely Replacement (MMLR) [22]. In [16], a Gabor dictionary is used in the reconstruction of thickness with a Matching Pursuit (MP)-based algorithm that penalizes a relaxed support measure corresponding to the p -pseudonorm with 0 < p < 1 .

3. Off-Grid Events and Dictionary Expansion

Aside from poor matrix conditioning, another problem known as off-grid deviation [23] limits the applicability of inverse-problem-based approaches on signal and image reconstruction. Figure 2a illustrates a grid of n = 9 modelled positions, represented by gray dots. As three events (represented by black dots) are located on modelled positions, the corresponding data vector c can be synthesized according to the acquisition models of Equations (2) and (3). The same does not hold when an off-grid event (represented by a red dot) is added: attempts to reconstruct the locations and amplitudes for the corresponding events may fail, causing artifacts and degradation on the reconstructed image.
Some formulations have been proposed for off-grid signal reconstruction, mainly within the framework of Compressive Sensing. In [24], the acquisition model considers a perturbation matrix summed column-wise to the (here referred to as H ) regular discrete model matrix. The formulation is applied to direction-of-arrival (DOA) estimation using the derivatives of the columns of H with respect to the sampled parameters as perturbation matrix. In [25], an adaptation of the Orthogonal Matching Pursuit (OMP) algorithm is proposed where the columns of the model matrix are iteratively updated in order to accommodate variations in the parameters of the PSFs. The algorithm is applied to pulse-Doppler radar. In [26], the problem of continuous line spectral estimation is approached with an algorithm based on the atomic norm minimization, which is solved via semi-definite programming. Similarly to the 1 minimization, the atomic norm minimization promotes sparse solutions. In [27], the regression problem uses a Total Least Squares (TLS) penalization with sparsity constraints. The motivation is that the “errors-in-variables” assumption of the TLS regression might be able to capture the mismatch between the model matrix and the acquired data. The method is then applied to cognitive radio sensing and DOA estimation.
Our approach relies on the framework of dictionary expansion. Each column h n of the discrete model H of Equation (3) is replaced by K columns [ b 1 ( n ) , , b K ( n ) ] = B ( n ) R ( M × K ) so that a data vector c resulting from the acquisition of an event located in the neighborhood of an n-th modelled position can be approximated by some linear combination of B ( n ) , i.e., by B ( n ) x ( n ) , where x ( n ) R K . As a result, an arbitrarily acquired c might be approximated as
c n = 1 N B ( n ) x ( n ) .
In the two-dimensional case, the neighborhood of the n-position is the region within ( x n ± 0.5 Δ x , z n ± 0.5 Δ z ) . This is represented in Figure 2b, where the nine modelled locations give place to nine neighborhoods (local ROIs).
Two forms of approximation are proposed in [23] for one-dimensional linear time-invariant (LTI) problems. The first one is the Taylor approximation, which relies on the fact that small shifts on a waveform can be well approximated by its Taylor expansion, i.e., by linearly combining the original waveform and its time derivatives. In this case, the column b 1 ( n ) is identical to the original atom h n and the columns b k ( n ) for k > 1 correspond to its ( k 1 ) -th time derivatives. The second is the Polar approximation, which is motivated by the fact that the continuous manifold M of an LTI system lies over a hypersphere on the M-dimensional data space [23]. The PSFs of the neighborhood of each n-th modelled position are approximated by an arc of a circle and the the column h n is replaced by three normal vectors with the directions of the center ( b 1 ( n ) ) and the two trigonometric components ( b 2 ( n ) and b 3 ( n ) ) of the circle. While the Taylor approximation can be done for any order K, in the Polar case K always equals 3.
An extension of the Basis Pursuit (BP) formulation [28], referred to as Continuous Basis Pursuit (CBP), is proposed in [23] for the recovery of the expanded coefficients { x ( n ) } 1 n N . For the sake of conciseness, from this point on, we will represent sets { x ( n ) } 1 n N simply as { x ( n ) } . The formulation of CBP is given by
{ x ^ ( n ) } = arg min { x ( n ) } 1 2 σ 2 c n = 1 N B ( n ) x ( n ) 2 2 + λ n = 1 N | x 1 ( n ) | ,
s . t . { x ( n ) } C ,
where the constraint set C prevents recovered expanded coefficients from having any arbitrary values that do not represent actual PSFs. The definition of the convex set C varies according to the type of approximation used. The 1 norm of a vector composed by the first element x 1 ( n ) of each K-tuple x ( n ) is used to obtain sparse solutions.
In [29], a low-rank approximation of the PSFs within the neighborhood of each n-th modelled position is performed by means of a Singular Value Decomposition (SVD). The continuous manifold drawn by τ in a local ROI is sampled with a very fine grid of R locations, generating R columns that form a matrix M ( n ) R M × R , as represented in Figure 2f. Each matrix M ( n ) then undergoes an SVD decomposition and the K first left singular vectors compose the corresponding expanded coefficients B ( n ) for the n-th local ROI.
An adaptation of the OMP [30] algorithm, referred to as Continuous OMP (COMP), is also presented in [29]. It aims at solving the 2 0 problem
{ x ^ ( n ) } = arg min { x ( n ) } ( x 1 ( 1 ) , , x 1 ( N ) ) 0 ,
s . t . c n = 1 N B ( n ) x ( n ) 2 2 ϵ { x ( n ) } C ,
where the symbol · 0 denotes the 0 pseudonorm, i.e., the cardinality (number of nonzero elements) of a vector.
In [31], a minimize–maximum (Minimax) formulation is presented for the definition of the expanded set { B ( n ) } . The resulting approximation minimizes the maximum residual within the representation of each n-th local ROI. It is motivated by the assumption that the off-grid deviation from a discrete grid follows a uniform distribution; therefore, the off-grid error should be as constant as possible, not privileging any distance from originally modelled positions.

4. Rank-K Approximation of Local Manifolds

Our criterion to determine B ( n ) is based on the SVD expansion, which has been proposed for one-dimensional, shift-invariant problems [29]. The extension to D-dimensional problems relies mainly on the first step of the process, which is a fine sampling of each local manifold M n : here, the regular, fine grid is defined for all D dimensions. This extension is facilitated by the fact that the formulation is non-parametric, i.e., the deviation from originally modelled positions is not mapped onto any independent variable and does not play any role on the definition on the bases. On the other hand, in the Taylor, Polar [23] and Minimax [31] expansions, the off-grid deviation is a parameter from which the elements of the expanded dictionary are derived. Consequently, except for the Taylor expansion, their extensions to two or higher dimensions are not promptly defined.

4.1. Highly Coherent Discrete Local Manifolds

Figure 2d shows an illustrative example of a D-manifold embedded in an M-dimensional data space. In this case, D = 2 and M = 3 . The nine D-dimensional modelled positions shown in Figure 2a correspond here to nine samples of the M-dimensional manifold, as well represented by gray dots in Figure 2d. The red dot corresponds to the data caused by the off-grid reflector from Figure 2a.
Figure 2e shows the same manifold as Figure 2d, but, instead of having N modelled positions, it divides the manifold into N local manifolds
M n : = { y ( x , z ) : x [ x n 0.5 Δ x , x n + 0.5 Δ x ] , z [ z n 0.5 Δ z , z n + 0.5 Δ z ] } ,
which correspond to the N local ROIs of Figure 2b.
We start by performing a fine sampling on each local manifold M n , as represented in Figure 2f. In practice, this means acquiring the PSF of a set of points from a fine grid of R points defined for each local ROI (Figure 2c). The result is a matrix M ( n ) R M × R , whose columns are local manifold samples. The finer this grid is, the better the continuous local manifold is represented by the discrete dataset M ( n ) . For simplicity of notation, we keep regular spacing δ x and δ z for the lateral and axial directions, respectively. The number of sampled points is R = R x × R z , where R x and R z are the number of locations defined on the lateral and axial directions, respectively. In the example of Figure 2c, R x = R z = 7 , thus R = 49 .
Our sampling includes the boundaries of the local ROIs. For this reason, the relation between the spacing and the number of locations on the lateral direction is given by
δ x = Δ x R x 1
and the same holds for the axial direction.

4.2. SVD Expansion

For each matrix M ( n ) , a rank-K approximation M ˜ ( n ) R M × R is to be defined and also factorized in the form
M ˜ ( n ) = B ( n ) F ( n ) ,
where B ( n ) is an orthonormal basis matrix and F ( n ) R K × R modulates B ( n ) to form M ˜ ( n ) . Any approximation creates a residual matrix R ( n ) R M × R defined by the difference
R ( n ) = M ( n ) B ( n ) F ( n ) .
The SVD expansion is defined by the minimization of the Frobenius norm [19] of R ( n ) :
B ^ ( n ) , F ^ ( n ) = arg min B ( n ) , F ( n ) M ( n ) B ( n ) F ( n ) F .
According to the Eckart–Young theorem, a solution for Equation (12) is achieved by a truncated SVD [32]. Consider the SVD of M
M ( n ) = U Σ V T ,
where U R M × R is the unitary matrix of left singular vectors, Σ R R × R is the diagonal matrix of singular values and V R N × R is the unitary matrix of right singular vectors [19]. The rank-K SVD truncation is obtained by using only the K largest singular values from Σ and the K corresponding vectors from U and V . This low rank approximation is given by
M ˜ ( n ) = U ˜ Σ ˜ V ˜ T ,
where U ˜ R M × K , Σ ˜ R K × K and V ˜ R R × K .
The K columns of U ˜ form an orthonormal basis for M ˜ ( n ) and composes the expanded set B ( n ) , while the product Σ ˜ V ˜ T compose the modulating matrix F ( n ) :
B ( n ) = U ˜ ,
F ( n ) = Σ ˜ V ˜ T .
Naturally, large values for K mean more degrees of freedom in the approximation, which reduces the residuals. Figure 3a shows how the value of K affects the Frobenius norm of R ( n ) for the centermost local ROI of the acquisition set presented in Section 6.1. The values of the 35 first singular values σ k are shown in Figure 3b. The 75 individual residual norms r i for K = 5 , 10 and 20 are shown in Figure 3c.
It shall be noted that the processes presented from Equation (12) to Equation (15b) have to be independently performed for every n-th local ROI. Although the construction of expanded dictionaries is computationally demanding, it is an offline procedure that is carried only once for each given acquisition set.

5. Reconstruction Algorithm

5.1. Limitations of Conic Constraints

Two main algorithms were proposed to work with expanded dictionaries: the convex CBP [23] and the greedy COMP [29]. The first one aims at solving the problem of Equation (6) while the second attempts to solve the problem of Equation (7). A hybrid approach called Interpolating Band-excluded Orthogonal Matching Pursuit (IBOMP) was also proposed and applied to frequency estimation (FE) and time delay estimation (TDE) [33]. Basically, it performs a rough greedy estimation of the support of the solution, followed by a refining convex optimization.
In order to implement a constraint set C , all the aforementioned algorithms have at least one step involving a constrained convex optimization where the constraints define either first-order (SVD, Minimax and Taylor) or second-order (Polar) cones. Figure 4a illustrates an example of a first-order cone for K = 2 . The black curved line represents the projection onto the basis B ( n ) of a continuous one-dimensional PSF manifold. The R vectors that compose a local manifold matrix M ( n ) , when projected onto B ( n ) , result in vectors f ( n ) , represented by the dots, which compose the columns of F ( n ) . When a reconstruction is performed, the recovered coefficients set x ( n ) R 2 for this n-th local ROI is constrained to lie within a first-order cone, represented by the shadowed area (which extends indefinitely to the right). This cone is defined by two linear constraints that impose an upper and a lower bound for the relation x 2 ( n ) / x 1 ( n ) , combined with a non-negativity constraint for the first component x 1 ( n ) . This constraint set aims to avoid arbitrary combinations for x ( n ) that do not represent positively-weighted copies of actual manifold samples. The upper black dot defines the upper angle of the cone, and is defined by the modulating matrix F ( n ) as max i ( f 2 , i ( n ) / f 1 , i ( n ) ) , i.e., the maximum relation between the first and second components found among the projections of M ( n ) . Similarly, the lower black dot is defined by min i ( f 2 , i ( n ) / f 1 , i ( n ) ) , and defines the lower angle of the cone. For higher orders of K, such a cone is defined for all K 1 relations between each k-th ( k 2 ) component and the first one. The resulting linear constraint set is defined as [29,31]
min 1 i R f k , i ( n ) f 1 , i ( n ) x k ( n ) x 1 ( n ) max 1 i R f k , i ( n ) f 1 , i ( n ) ,
f 1 , i ( n ) 0 ,
k { 2 , , K } , n { 1 , 2 , , N } ,
where f k , i ( n ) denotes the element on the k-th line and i-th column on F ( n ) . The principle is similar for the Polar expansion, though in that case the cones are of second order [23].
Notice that the cone-based convex constraints assume that the projection of M ( n ) on the K components of B ( n ) yields relatively large, positive, small-variance values for the first component and small values for the remaining, yielding relatively small values for minimum and maximum relations of Equation (16c). If this assumption is broken, the cone will span too large an area of the right half-plane, i.e., it will constrain less, being less effective, as represented in Figure 4b. In some cases, defining the the cone is not even possible, as depicted in Figure 4c.
Assuring a well behaved relation between the first and the remaining components, as shown in Figure 4a, implies choosing considerably small values for Δ x and Δ z , which limits the applicability of recovery algorithms based on conic constraints. For instance, on the simulated acquisition set of Section 6.1, choosing Δ x = Δ z = 0.2 mm still causes the first component to have both positive and negative values on certain local manifolds, which makes the CBP [23], COMP [29] and IBOMP [33] algorithms not applicable.

5.2. Non-Convex Constraints

Instead of using conic constraints, our algorithm attempts to constrain each K-tuple of recovered coefficients x ( n ) to be similar to any column of the modulating matrix F ( n ) . We translate “similarity” as high correlation, as formalized in the non-convex constraint set
max 1 i R x ( n ) , f i ( n ) x ( n ) f i ( n ) μ c , n { 1 , 2 , , N } ,
where a , b = a T b denotes the inner product of two vectors.
The minimum correlation parameter μ c controls how similar to any of the manifold samples on M ( n ) a recovered event must be. If a given x ( n ) passes the test of Equation (17), proving to be sufficiently similar to some i-th modulating vector f i ( n ) , then the approximation
m ˜ i ( n ) x ( n ) f i ( n ) = B ( n ) f i ( n ) x ( n ) f i ( n ) B ( n ) x ( n )
is assumed and the product B ( n ) x ( n ) is considered as a valid weighted copy of a PSF within the n-th local ROI, rather than an arbitrary combination of the n-th basis vectors. This constraint is imposed by our greedy algorithm on the decision of which expanded set B ( n ) will be added to the reconstruction problem at each iteration.

5.3. OMP for Expanded Dictionaries

The proposed algorithm, summarized in Algorithm 1, is an extension of the OMP algorithm, referred to as OMP for Expanded Dictionaries (OMPED). It attempts to solve a problem similar to Equation (7) with the non-convex constraint set C defined in Equation (17). The stop criterion is based on the residual yielded by the LS solution with a given cardinality, yet, instead of comparing the residual to a fixed parameter ϵ , we compare it to an estimate of the current residual that takes into account the expected acquisition noise and the estimated residuals resulting from the reduced-rank approximation.
Algorithm 1 OMP for Expanded Dictionaries (OMPED)
Input: { B ( n ) } , { F ( n ) } , { R ( n ) } , c , e noise , μ c , Δ μ
1:
S
2:
e c
3:
repeat
4:
j Compute from Equation (21)
5:
while j = do
6:
   μ c μ c Δ μ
7:
   j Compute from Equation (21)
8:
end while
9:
S S { j }
10:
{ x ( n ) } Compute from Equation (22b)
11:
e Compute from Equation (23)
12:
e rank Compute from Equation (24)
13:
e est Compute from Equation (25)
14:
until e est e 2 or S C =
Output:S, { x ( n ) } n S
The input parameter e noise contains the expected 2 norm of the acquisition noise. In practice, this value can be obtained from acquisitions with samples of the inspected material known to have neither discontinuities nor other sort of scatterers. For our simulations, we use the relation
e noise 2 = w 2 2 M σ 2 ,
which holds if the noise vector w contains white Gaussian noise with variance σ 2 . The approximation of Equation (19) becomes an equality as M . We assume the equality and use e noise = M σ 2 .
We define the support S of the solution, which is initialized as the empty set, and its complement S c = { 1 , , N } S . The solution residual e R M is initialized with the vector of acquired data c on line 2.
At each iteration, an index j S c is added to S as we choose the expanded set B ( j ) that is capable of causing the maximal decrease on the energy of the residual, as represented on the left side of Equation (20). Since the columns of each B ( n ) are orthonormal, the identity
j ^ = arg min j e B ( j ) B ( j ) T e 2 = arg max j B ( j ) T e 2
holds as a consequence of Parseval’s relation [34], which allows us to perform the simpler operation of taking the norm of each product B ( j ) T e .
This operation is a generalization of the measurement of maximum correlation on the original OMP [30]. A constraint based on Equation (17) is imposed to prune candidates that do not accomplish the minimum correlation required. The resulting criterion is formalized as
j ^ = arg max j S C B ( j ) T e 2 s . t . max 1 i R B ( j ) T e , f i ( j ) B ( j ) T e f i ( j ) μ c .
The constraint in Equation (21) allows for the recovery of only positive-amplitude events. It can be adapted to consider both positive and negative amplitudes by simply replacing the inner product by its absolute value | B ( j ) T e , f i ( j ) | .
The algorithm considers the case where no index meets the correlation criterion of Equation (21). This case is treated from line 5 to line 8: while the problem of Equation (11) remains infeasible, a decrease of Δ μ is made on the parameter μ c and a new attempt to compute the index j is performed.
The support S is then updated to include the new index j (line 9) and is used to compute the coefficients
{ x ^ ( n ) } = arg min { x ( n ) } c n = 1 N B ( n ) x ( n ) 2 2 ,
s . t . x ( n ) = 0 , n S c ,
where 0 R K is the zero vector. The updated residual is yielded as
e = c n S B ( n ) x ( n ) .
Were the manifold approximation exact, e in Equation (23) would be composed strictly of: (1) PSFs located at local ROIs with the corresponding indices not yet added to the support S and (2) additive noise. In that case, we could use the widespread stop criterion that compares e 2 to the expected noise power. However, our residual estimate must take into account the rank-K approximation. This estimate is computed on vector e rank R M as
e rank = n S r i ^ ( n ) x ( n ) f i ^ ( n ) ,
where i ^ = arg max 1 i R x ( n ) , f i ( n ) x ( n ) f i ( n ) ,
and r i ( n ) denotes the i-th column from R ( n ) . Based on Equation (17), the index i in Equation (24b) is a function of n: for every index n in the current support S, the correlations performed in (24b) estimate which i-th PSF within the n-th local manifold best explains the recovered coefficients x ( n ) (see Figure 2c,f). The residual r i ( n ) , from the dictionary low-rank approximation, is then used as a template for the estimation of the current approximation residual. The amplitude estimate is taken from the ratio between the norms of the recovered coefficients x ( n ) and of the similar modulating vector f i ( n ) .
The current total residual norm is estimated as
e est = ( e rank 2 2 + e noise 2 ) 1 2 ,
where the summation is performed under the assumption that the acquisition noise and the vector e rank have negligible correlation.
The algorithm greedily increases the support until the estimated residual norm e est reaches the norm e of the actual residual yielded by the LS or all indices n = 1 , , N have been added to the support S.

5.4. Recovery of Locations and Amplitudes

Recalling the approximation m i ( n ) B ( n ) f i ( n ) , we determine i by finding out which f i ( n ) most correlates to x ( n ) :
i ^ ( n ) = arg max 1 i R x ( n ) , f i ( n ) x ( n ) f i ( n ) , n S .
The amplitude estimations v n result form the ratios between the norms of x ( n ) and of the chosen template f i ( n ) :
v n = x ( n ) f i ( n ) , n S , i as in Equation (26).
As consequence, the spatial resolution of the reconstructed events equals the fine sampling represented in Figure 2c, i.e., δ x and δ z for the lateral and axial axes, respectively.

6. Empirical Results

6.1. Simulated Acquisition Set

To simulate the ultrasound NDT acquisition set from [21], represented in Figure 5a, we used the Field II package [15] for Matlab® (The MathWorks Inc., Natick, MA, USA). A piston transducer with 3 mm radius (125 μ m mathematical element size) interrogates a steel sample object (sound speed c = 5680 m/s). The excitation pulse has center frequency f c = 5 MHz and 6 dB fractional bandwidth of 100 % . The simulated transducer slides horizontally along the surface of the object, acquiring scanlines from 31 lateral positions u i , from u 0 = 0 mm to u 30 = 30 mm (center of transducer), with a distance of 1 mm between consecutive lateral positions. The 31 scanlines are sampled with sampling rate f s = 25 MHz and concatenated to form the acquisition vector c .
Following [21], the model grid has 31 × 41 = 1271 modelled locations distributed with regular spacing of 1 mm on both x- and z-directions. On the x-direction, the locations are the same as the transducer positions, i.e., x = 0 , 1 mm , , 30 mm . On the z-direction, 41 locations are modelled regularly between 18 mm and 58 mm, i.e., z = 18 mm , 19 mm , , 58 mm .
As explained in Section 4.1, in the expanded acquisition model, the grid locations give place to local ROIs. Our expanded model has 1271 local ROIs with Δ x = Δ z = 1 mm, with centers corresponding to the modelled locations of the regular model. Consequently, our ROI extends from x = 0.5 mm to x = 30.5 mm and from z = 17.5 mm to z = 58.5 mm. The highly coherent local manifolds were created with R x = 5 and R z = 15 , thus R = 75 . Therefore, δ x = 250 μ m and δ z = 71.4 μ m.
We simulated the acquisition for 200 cases of five unity amplitude scatterers randomly distributed over the ROI. The scatterers’ positions were not forced over any kind of grid. White Gaussian noise with three different levels ( σ = 0 , 0.08 , 0.12 ) was added to each simulated acquisition. Since the energy of the acquired signal (without noise) varies according to factors such as distance to transducer and constructive/destructive interference, we consider that the parametrization of noise in terms of its standard deviation σ is more appropriate than signal-to-noise ratio (SNR). To provide a visual notion of the noise levels, Figure 5b shows an extract of acquired data for the three noise levels from an acquisition where a single scatterer was placed on the center of the ROI. Scanlines from the three centermost positions of the tranducer are concatenated.

6.2. Recovery Accuracy

To compute the accuracy on the recovery of scatterers, we ran OMPED with a fixed number of five iterations, with μ c = 0.8 , Δ μ = 0.1 and K varying from 2 to 10 for the 200 simulated acquisitions with the three levels of noise. Each recovered scatterer distant less than 0.5 mm in both axial and lateral directions from the closest original simulated scatterer was computed as a hit—otherwise, it was computed as a miss. Figure 6a shows the percentage of misses from 1000 recovered scatterers for all nine values of K and three noise levels. Even for the highest level of noise, misses kept below 10 % for 6 K 10 . For comparison, we ran OMP with the regular dictionary H on the same set of simulated acquisitions. The resulting percentages of misses were 38.9 % , 42.4 % and 45.2 % for the noise levels σ = 0 , 0.08 and 0.12 , respectively.
A small increase in the count of misses is observed for values of K 8 . This is possibly explained by the fact that, for K 8 , increasing K adds few useful information to the dictionary at the cost of increasing coherence. For the SVD basis, the value of the singular values σ k can be used as a measure of useful information. Figure 3b shows how σ k behaves for the centermost local manifold M ( 636 ) . Notice that values of σ k for k 8 are significantly smaller than the previous ones.
For every hit, the distance between the original and the recovered scatterers was computed. The average distances are shown in Figure 6b.
The computation of hits and misses does not take into account the amplitude of recovered scatterers, i.e., recovered scatterers are implicitly considered as having unity amplitude. To endorse this assumption, the average amplitudes of recovered events are shown in Figure 6c, where the bars indicate one standard deviation above and below the average. Notice that, for all cases, the average amplitudes are between 0.98 and 1.01 , i.e., the average amplitude error is less than 2 % . The average absolute amplitude resulting from the reconstructions with OMP using the regular dictionary H were 0.70 , 0.70 and 0.71 for the noise levels σ = 0 , 0.08 and 0.12 , respectively.

6.3. Estimation of Residual and Stop Criterion

To assess the accuracy of the stop criterion, OMPED was executed one more time on the 5-scatterer dataset of Section 6.1, this time with the residual-based stop criterion defined on line 14 of Algorithm 1, with a maximum of 10 iterations. Because all images contained five scatterers, the algorithm was expected to stop at the 5th iteration. The histogram of Figure 7a shows this outcome: the peak of occurrences is on iteration 5. The frequencies on the neighboring final iterations 4 and 6 are also sensibly greater than on the remaining iterations (except for the maximum 10). The maximum iteration allowed was 10, at which the algorithm stopped when e est failed to reach e . The results for values of K from 2 to 10 are summed on the histogram of Figure 7a. A total of 5400 reconstruction (3 noise levels × 200 images × 9 orders K) are computed.
Figure 7b shows an example of the evolution of the regression residual norm e and the estimated residual norm e est . As new events are iteratively added to the solution, the latter decreases while the former increases. On iteration 5, e drops below e est and OMPED correctly meets the stop criterion, yielding a final solution with cardinality 5. White Gaussian noise with σ = 0.12 was added to the data. OMPED was ran with SVD ( K = 8 ) dictionary.

6.4. Reconstructed Images: Examples

Figure 8a shows the ground truth for a simulation from the dataset of Section 6.1. Gaussian noise was added to the acquired data with σ = 0.08 . The reconstructed image using OMPED with SVD dictionary ( K = 8 ) is shown in Figure 8b. No limit was imposed on the number of iterations, i.e., the algorithm correctly stopped at the 5th iteration based on the values of the estimated and actual residuals. The activated pixels are the same on the ground truth of Figure 8a and on the OMPED result of Figure 8b. While all simulated scatterers had unity amplitude, the recovered amplitudes ranged from 0.9398 to 1.0387 . Both Figure 8a,b have 41 × 31 pixels corresponding to the local ROIs of the expanded model.
The result of the reconstruction using OMP with the regular dictionary model H is shown in Figure 8c. We ran seven iterations of the algorithm in order to show that, beyond iteration 4, the algorithm created artifacts around the leftmost scatterer instead of identifying the bottom-right scatterer. The recovered amplitudes also display small and even negative values (the image shows absolute, normalized values). Moreover, the bottom-left scatterer is displaced one pixel to the left on the reconstructed image.
Figure 8d shows the image yielded by the LS (unregularized) solution of Equation (4). As is common in unregularized model-based solutions, the image is dominated by noise [35]. We also applied 1 regularization to the LS problem, which corresponds to the BP formulation [28]
v ^ = arg min v c H v 2 2 + λ v 1 .
The 1 -regularized formulation was solved with L1_LS package for Matlab [36]. The resulting image is shown in Figure 8e. While a small value for λ yields an image dominated by noise, such as that of Figure 8d, larger values cause the image to be too sparse, suppressing some features. This is a consequence of the penalization of recovered amplitudes on Equation (28). The chosen regularization parameter λ = 2.0691 minimizes the norm v v ^ 2 , where v is the ground truth and v ^ is the BP result.

7. Conclusions

To cope with the problem of off-grid deviation in image reconstruction from pulse-echo ultrasound data, we developed a technique of dictionary expansion based on a highly coherent sampling of the PSF manifold followed by a rank reduction procedure, as well as a generalization of the OMP algorithm with non-convex constraints. Based on [29], the criterion for the rank reduction is the minimization of the Frobenius norm of the resulting residuals.
Since no assumption is made regarding the geometry of the continuous PSF manifold, our expansion formulation is applicable to both shift-invariant and shift-variant problems. On the other hand, for instance, the Polar expansion [23] is conceived based on the fact that the PSF manifold of any shift-invariant system lies over a hypersphere. In a two-dimensional ultrasound (our main motivating application), the fact that the Spatial Impulse Response (SIR) is spatially variant [15,37] puts the direct acquisition model in the class of shift-variant systems.
The criterion for definition of the order K of expansion may vary according to each application. In cases where it is possible to carry out simulations (as presented here) or a relevant amount of data with accessible ground truth is available, K can be determined empirically. Moreover, in our case, a minimum in the number of misses is identifiable and lies close to a transition on the baseline of singular values shown in Figure 3b. A suggestion for future studies is the development of a generalized criterion for the definition of K. The behavior of the singular values of matrices M ( n ) is potentially a starting point for such investigation.
The original OMP algorithm [30] is a particular case of OMPED where K = 1 and the parameter μ c (Equations (17) and (21)) is set to an arbitrarily large negative value. In both OMP and OMPED, the residual vector e on each iteration is orthogonal to all active elements of the dictionary, which places OMPED in the family of Orthogonal Matching Pursuit algorithms. The same does not hold for the COMP algorithm presented in [29]: the fact that the LS regression performed at each iteration contains linear constraints may result in eventual coherence between the residual and the active atoms.
Another particularity of OMPED in regard to previously proposed algorithms for expanded dictionaries [23,29,33] is that it is not based on conic constraints, which removes any restrictions on the choice of the sizes Δ x and Δ z (and further dimensions if that is the case) for the division of the ROI into local ROIs.
The adaptation of OMP into OMPED, with a constraint imposed on the selection of the index added the support at each iteration, might be replicable to other greedy search algorithms. The class of forward-backward algorithms is of special interest in signal and image recovery because of its capacity of later “correction” of “wrong” choices made on the selection of indices to add to the support [38,39], which constitutes a motivation for future investigation.
The computation of the estimated residual e est on OMPED may be subject to improvement in order to increase the accuracy of the stop criterion (see Figure 7a). Decreasing the variance of the residuals r i ( n ) caused by the low-rank approximation inside each local ROI (i.e., flattening the surfaces of Figure 3c) would cause the inaccuracies on the computation of high resolution locations to have a smaller impact on the computation of e est . This may be achieved with a different criterion for the rank reduction than the LS. For instance, an extension of the Minimax dictionary expansion [31].
One limitation of our technique is that one single point-like event is identifiable inside each local ROI. The search for a means to overcome this limitation, allowing for the recovery of several scatterers inside the same local ROI, is a relevant topic for further investigation and may broaden the applicability of the proposed technique.
Finally, our simulated data considered point-like reflectors, with spatial coordinates ( x , z ) as the only nonlinear parameters. The ultrasound NDT literature contains parametric reflection models for more complex discontinuity structures, such as spherical voids and circular cracks, where the distortion of ultrasound waves is modelled as a nonlinear function of parameters like diameter and angle to the surface [40,41]. The proposed method is applicable to those cases as long as those parameters are comprised in the parameter set τ in Equation (1) and sampled like the parameters of spatial location. In this case, characterization of discontinuities could be performed along with location. Classification of discontinuities could also be jointly performed if dictionaries for several types of discontinuities are combined. An equivalent principle has been used in the joint detection and identification of neuron activity using SVD [29] and Taylor [42] expanded dictionaries.

Author Contributions

All authors are part of the Lineup for Inverse Problems, Signal and Image Processing and Reconstruction (LIPRO) at the Graduate Program in Electrical and Computer Engineering (CPGEI), at the Federal University of Technology-Paraná (UTFPR), Brazil. T.A.R.P. was responsible for generating the datasets, writing and running the MATLAB codes, and describing the models and methods. M.V.W.Z. and D.R.P. were responsible for reviewing the mathematical accuracy of methods, organizing the paper argument and reviewing the ultrasonic model accuracy.

Funding

The authors would like to thank Petrobras/CENPES and CNPq grant 312023/2015-4 for finacial support.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Fessler, J.A. Model-based image reconstruction for MRI. IEEE Signal Process. Mag. 2010, 27, 81–89. [Google Scholar] [CrossRef] [PubMed]
  2. Ollinger, J.M.; Fessler, J.A. Positron-emission tomography. IEEE Signal Process. Mag. 1997, 14, 43–55. [Google Scholar] [CrossRef] [Green Version]
  3. Censor, Y. Finite series-expansion reconstruction methods. Proc. IEEE 1983, 71, 409–419. [Google Scholar] [CrossRef]
  4. Lavarello, R.; Kamalabadi, F.; O’Brien, W.D. A regularized inverse approach to ultrasonic pulse-echo imaging. IEEE Trans. Med. Imaging 2006, 25, 712–722. [Google Scholar] [CrossRef] [PubMed]
  5. Zanin, L.; Zibetti, M.; Schneider, F. Conjugate gradient and regularized Inverse Problem-Based solutions applied to ultrasound image reconstruction. In Proceedings of the 2011 IEEE International Ultrasonics Symposium (IUS), Orlando, FL, USA, 18–21 October 2011; pp. 377–380. [Google Scholar]
  6. Zanin, L.G.S.; Schneider, K.F.; Zibetti, M.V.W. Regularized Reconstruction of Ultrasonic Imaging and the Regularization Parameter Choice. In Proceedings of the International Conference on Bio-Inspired Systems and Signal Processing, Vilamoura, Portugal, 1–4 February 2012; pp. 438–442. [Google Scholar]
  7. Desoky, H.; Youssef, A.B.M.; Kadah, Y.M. Reconstruction using optimal spatially variant kernel for B-mode ultrasound imaging. In Proceedings of the SPIE Medical Imaging 2003: Ultrasonic Imaging and Signal Processing, San Diego, CA, USA, 23–28 February 2003; pp. 147–153. [Google Scholar]
  8. Lingvall, F.; Olofsson, T.; Stepinski, T. Synthetic aperture imaging using sources with finite aperture: Deconvolution of the spatial impulse response. J. Acoust. Soc. Am. 2003, 114, 225–234. [Google Scholar] [CrossRef] [PubMed]
  9. Lingvall, F.; Olofsson, T. On time-domain model-based ultrasonic array imaging. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 2007, 54, 1623–1633. [Google Scholar] [CrossRef] [PubMed]
  10. Olofsson, T.; Wennerstrom, E. Sparse Deconvolution of B-Scan Images. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 2007, 54, 1634–1641. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  11. Viola, F.; Ellis, M.A.; Walker, W.F. Time-domain optimized near-field estimator for ultrasound imaging: Initial development and results. IEEE Trans. Med. Imaging 2008, 27, 99–110. [Google Scholar] [CrossRef] [PubMed]
  12. Schmerr, L.W., Jr. Fundamentals of Ultrasonic Phased Arrays; Springer: Berlin, Germany, 2014; Volume 215. [Google Scholar]
  13. Smith, N.; Webb, A. Introduction to Medical Imaging: Physics, Engineering and Clinical Applications; Cambridge Texts in Biomedical Engineering, Cambridge University Press: Cambridge, UK, 2010. [Google Scholar]
  14. Jensen, J.A. A model for the propagation and scattering of ultrasound in tissue. Acoust. Soc. Am. 1991, 89, 182–190. [Google Scholar] [CrossRef] [Green Version]
  15. Jensen, J. Simulation of advanced ultrasound systems using Field II. In Proceedings of the 2004 IEEE International Symposium on Biomedical Imaging: Nano to Macro, Arlington, VA, USA, 15–18 April 2004; Volume 1, pp. 636–639. [Google Scholar] [CrossRef]
  16. Mor, E.; Azoulay, A.; Aladjem, M. A matching pursuit method for approximating overlapping ultrasonic echoes. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 2010, 57, 1996–2004. [Google Scholar] [CrossRef] [PubMed]
  17. Carcreff, E.; Bourguignon, S.; Idier, J.; Simon, L. A linear model approach for ultrasonic inverse problems with attenuation and dispersion. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 2014, 61, 1191–1203. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  18. Kruizinga, P.; van der Meulen, P.; Fedjajevs, A.; Mastik, F.; Springeling, G.; de Jong, N.; Bosch, J.G.; Leus, G. Compressive 3D ultrasound imaging using a single sensor. Sci. Adv. 2017, 3. [Google Scholar] [CrossRef] [PubMed]
  19. Golub, G.H.; Van Loan, C.F. Matrix Computations (Johns Hopkins Studies in Mathematical Sciences), 3rd ed.; The Johns Hopkins University Press: Baltimore, MD, USA, 1996. [Google Scholar]
  20. Hansen, C. Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion; SIAM: Philadelphia, PA, USA, 1998. [Google Scholar]
  21. Guarneri, G.A.; Pipa, D.R.; Junior, F.N.; de Arruda, L.V.R.; Zibetti, M.V.W. A Sparse Reconstruction Algorithm for Ultrasonic Images in Nondestructive Testing. Sensors 2015, 15, 9324–9343. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. Chi, C.Y.; Goutsias, J.; Mendel, J. A fast maximum-likelihood estimation and detection algorithm for Bernoulli-Gaussian processes. In Proceedings of the ICASSP ’85 IEEE International Conference on Acoustics, Speech, and Signal Processing, Tampa, FL, USA, 26–29 April 1985; Volume 10, pp. 1297–1300. [Google Scholar] [CrossRef]
  23. Ekanadham, C.; Tranchina, D.; Simoncelli, E.P. Recovery of Sparse Translation-Invariant Signals With Continuous Basis Pursuit. IEEE Trans. Signal Process. 2011, 59, 4735–4744. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Yang, Z.; Zhang, C.; Xie, L. Robustly stable signal recovery in compressed sensing with structured matrix perturbation. IEEE Trans. Signal Process. 2012, 60, 4658–4671. [Google Scholar] [CrossRef]
  25. Teke, O.; Gurbuz, A.C.; Arikan, O. A robust compressive sensing based technique for reconstruction of sparse radar scenes. Dig. Signal Process. 2014, 27, 23–32. [Google Scholar] [CrossRef] [Green Version]
  26. Tang, G.; Bhaskar, B.N.; Shah, P.; Recht, B. Compressed Sensing Off the Grid. IEEE Trans. Inf. Theory 2013, 59, 7465–7490. [Google Scholar] [CrossRef] [Green Version]
  27. Zhu, H.; Leus, G.; Giannakis, G.B. Sparsity-Cognizant Total Least-Squares for Perturbed Compressive Sampling. IEEE Trans. Signal Process. 2011, 59, 2002–2016. [Google Scholar] [CrossRef] [Green Version]
  28. Chen, S.S.; Donoho, D.L.; Saunders, M.A. Atomic Decomposition by Basis Pursuit. SIAM J. Sci. Comput. 1998, 20, 33–61. [Google Scholar] [CrossRef] [Green Version]
  29. Knudson, K.C.; Yates, J.; Huk, A.; Pillow, J.W. Inferring sparse representations of continuous signals with continuous orthogonal matching pursuit. Adv. Neural Inf. Process. Syst. 2014, 27, 1215–1223. [Google Scholar]
  30. Tropp, J.A.; Gilbert, A.C. Signal Recovery From Random Measurements Via Orthogonal Matching Pursuit. IEEE Trans. Inf. Theory 2007, 53, 4655–4666. [Google Scholar] [CrossRef] [Green Version]
  31. Passarin, T.A.R.; Pipa, D.R.; Zibetti, M.V.W. A minimax dictionary expansion for sparse continuous reconstruction. In Proceedings of the 2017 25th European Signal Processing Conference (EUSIPCO), Kos, Greece, 28 August–2 September 2017; pp. 2136–2140. [Google Scholar] [CrossRef]
  32. Eckart, C.; Young, G. The approximation of one matrix by another of lower rank. Psychometrika 1936, 1, 211–218. [Google Scholar] [CrossRef]
  33. Fyhn, K.; Marco, F.; Holdt, S. Compressive parameter estimation for sparse translation-invariant signals using polar interpolation. IEEE Trans. Signal Process. 2015, 63, 870–881. [Google Scholar] [CrossRef]
  34. Barrett, H.; Myers, K. Foundations of Image Science; Wiley Series in Pure and Applied Optics; Wiley-Interscience: Hoboken, NJ, USA, 2004. [Google Scholar]
  35. Bovik, A. Handbook of Image & Video Processing; Academic Press Series in Communications, Networking and Multimedia; Academic Press: Cambridge, MA, USA, 2000. [Google Scholar]
  36. Kim, S.J.; Koh, K.; Lustig, M.; Boyd, S.; Gorinevsky, D. An Interior-Point Method for Large-Scale L1-Regularized Least Squares. IEEE J. Sel. Top. Signal Process. 2007, 1, 606–617. [Google Scholar] [CrossRef]
  37. Tupholme, G.E. Generation of acoustic pulses by baffled plane pistons. Mathematika 1969, 16, 209–224. [Google Scholar] [CrossRef]
  38. Miller, A. Subset Selection in Regression; CRC Press: Boca Raton, FL, USA, 2002. [Google Scholar]
  39. Soussen, C.; Idier, J.; Brie, D.; Duan, J. From Bernoulli Gaussian Deconvolution to Sparse Signal Restoration. IEEE Trans. Signal Process. 2011, 59, 4572–4584. [Google Scholar] [CrossRef] [Green Version]
  40. Schmerr, L.W.; Song, S.J. Ultrasonic Nondestructive Evaluation Systems: Models and Measurements; Springer: Berlin, Germany, 2007. [Google Scholar]
  41. Velichko, A.; Bai, L.; Drinkwater, B. Ultrasonic defect characterization using parametric-manifold mapping. Proc. R. Soc. A 2017, 473, 20170056. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Ekanadham, C.; Tranchina, D.; Simoncelli, E.P. A unified framework and method for automatic neural spike identification. J. Neurosci. Methods 2014, 222, 47–55. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1. Acquisition of the point spread function (PSF). For each position ( x , z ) of the unity amplitude scatterer within the region of interest (ROI) (left side), an M-sample response y ( x , z ) R M (arranged as an M-pixel image on right side) is generated by the acquisition model. The set of all possible PSFs within the region of interest form a manifold M onto the data space.
Figure 1. Acquisition of the point spread function (PSF). For each position ( x , z ) of the unity amplitude scatterer within the region of interest (ROI) (left side), an M-sample response y ( x , z ) R M (arranged as an M-pixel image on right side) is generated by the acquisition model. The set of all possible PSFs within the region of interest form a manifold M onto the data space.
Sensors 18 04097 g001
Figure 2. (a) an illustrative discrete acquisition model with N = 3 × 3 = 9 modelled positions, represented by the gray dots. The black dots represent three well located events and the red dot represents an off-grid event. Because of the latter, the corresponding acquisition data vector c cannot be synthesized as a linear combination of the columns of the discrete model matrix H ; (b) the region of interest (ROI) is divided into N local ROIs with area Δ x × Δ z ; (c) each local ROI is sampled with a fine grid with lateral and axial distances δ x and δ z ; (d) on the space R M of acquired data, the set of all possible point spread function (PSFs) within the ROI form a manifold M . The gray dots are the PSFs of the modelled positions of Figure 2a. The black dots are on the grid, while the red dot is off-grid; (e) as the ROI is divided into N local ROIs (Figure 2b), the manifold is divided into N corresponding local manifolds; (f) the acquisitions over the fine grid on each n-th local ROI create R samples from the corresponding local manifold. Those samples compose matrix M ( n ) R M × R .
Figure 2. (a) an illustrative discrete acquisition model with N = 3 × 3 = 9 modelled positions, represented by the gray dots. The black dots represent three well located events and the red dot represents an off-grid event. Because of the latter, the corresponding acquisition data vector c cannot be synthesized as a linear combination of the columns of the discrete model matrix H ; (b) the region of interest (ROI) is divided into N local ROIs with area Δ x × Δ z ; (c) each local ROI is sampled with a fine grid with lateral and axial distances δ x and δ z ; (d) on the space R M of acquired data, the set of all possible point spread function (PSFs) within the ROI form a manifold M . The gray dots are the PSFs of the modelled positions of Figure 2a. The black dots are on the grid, while the red dot is off-grid; (e) as the ROI is divided into N local ROIs (Figure 2b), the manifold is divided into N corresponding local manifolds; (f) the acquisitions over the fine grid on each n-th local ROI create R samples from the corresponding local manifold. Those samples compose matrix M ( n ) R M × R .
Sensors 18 04097 g002
Figure 3. Approximation metrics for the centermost local ROI of the ultrasound acquisition set described in Section 6.1, with R = 75 ( R x = 5 and R z = 15 ). (a) Frobenius norm R ( n ) F of the residual matrix as a function of the order of approximation K; (b) 35 first singular values σ k from the singular value decomposition of M ( n ) ; (c) individual residual norms r i ( n ) 2 (of columns of R ( n ) ), spatially arranged according to the corresponding positions on the local ROI. The three surfaces correspond to K = 5 (top), K = 10 (middle) and K = 20 (bottom).
Figure 3. Approximation metrics for the centermost local ROI of the ultrasound acquisition set described in Section 6.1, with R = 75 ( R x = 5 and R z = 15 ). (a) Frobenius norm R ( n ) F of the residual matrix as a function of the order of approximation K; (b) 35 first singular values σ k from the singular value decomposition of M ( n ) ; (c) individual residual norms r i ( n ) 2 (of columns of R ( n ) ), spatially arranged according to the corresponding positions on the local ROI. The three surfaces correspond to K = 5 (top), K = 10 (middle) and K = 20 (bottom).
Sensors 18 04097 g003
Figure 4. (a) illustrative case of projection of local manifold samples M ( n ) on a basis B ( n ) , for K = 2 . The curved line represents the projection of a continuous one-dimensional manifold, while the dots represent the projection of the samples (columns of M ( n ) ) on B ( n ) . When Δ is sufficiently small, the projections have single-signed, relatively large values on the first component f 1 ( n ) and smaller values on the remaining components. In this case, the definition of a first-order cone (represented by the shadowed region) is possible and can be used in the reconstruction algorithm combined with a non-negativity constraint for the first component, ensuring that the recovered coefficients represent weighted copies of the local manifold, rather than other arbitrary combinations. The upper and lower angles of the cone depend on max i ( f 2 , i ( n ) / f 1 , i ( n ) ) and min i ( f 2 , i ( n ) / f 1 , i ( n ) ) , respectively; (b) as Δ increases, the angle of the cone may as well increase, making the constraint less effective, as a broader area is allowed for the recovered coefficients f ( n ) ; (c) an example where the definition of a convex cone is no longer possible. This imposes a limit on the definition of Δ .
Figure 4. (a) illustrative case of projection of local manifold samples M ( n ) on a basis B ( n ) , for K = 2 . The curved line represents the projection of a continuous one-dimensional manifold, while the dots represent the projection of the samples (columns of M ( n ) ) on B ( n ) . When Δ is sufficiently small, the projections have single-signed, relatively large values on the first component f 1 ( n ) and smaller values on the remaining components. In this case, the definition of a first-order cone (represented by the shadowed region) is possible and can be used in the reconstruction algorithm combined with a non-negativity constraint for the first component, ensuring that the recovered coefficients represent weighted copies of the local manifold, rather than other arbitrary combinations. The upper and lower angles of the cone depend on max i ( f 2 , i ( n ) / f 1 , i ( n ) ) and min i ( f 2 , i ( n ) / f 1 , i ( n ) ) , respectively; (b) as Δ increases, the angle of the cone may as well increase, making the constraint less effective, as a broader area is allowed for the recovered coefficients f ( n ) ; (c) an example where the definition of a convex cone is no longer possible. This imposes a limit on the definition of Δ .
Sensors 18 04097 g004
Figure 5. (a) simulated set (figure adapted from [21]). The transducer, fixed vertically at z = 0 , slides horizontally over the surface of the interrogated object, acquiring scanlines at 31 positions x = { u 0 , , u 30 } , corresponding to 0 mm up to 31 mm with 1 mm step. The scanlines are concatenated to form the acquired vector c . A PSF y ( x , z ) is determined by placing a unity amplitude scatterer on position ( x , z ) and acquiring the corresponding c ; (b) extracts from the acquired data for the three centermost transducer positions, with a unity amplitude scatterer located at the center of the ROI. White Gaussian noise was added with σ = 0 (up), σ = 0.08 (middle) and σ = 0.12 (bottom).
Figure 5. (a) simulated set (figure adapted from [21]). The transducer, fixed vertically at z = 0 , slides horizontally over the surface of the interrogated object, acquiring scanlines at 31 positions x = { u 0 , , u 30 } , corresponding to 0 mm up to 31 mm with 1 mm step. The scanlines are concatenated to form the acquired vector c . A PSF y ( x , z ) is determined by placing a unity amplitude scatterer on position ( x , z ) and acquiring the corresponding c ; (b) extracts from the acquired data for the three centermost transducer positions, with a unity amplitude scatterer located at the center of the ROI. White Gaussian noise was added with σ = 0 (up), σ = 0.08 (middle) and σ = 0.12 (bottom).
Sensors 18 04097 g005
Figure 6. (a) percentage of misses (from 1000 simulated events) as a function of K, for three levels of noise, with OMPED running with a fixed number of five iterations (each of the 200 simulated acquisition had five scatterers). Each recovered scatterer distant more than 0.5 mm in any direction (axial or lateral) from the closest original simulated scatterer was computed as a miss. A minimum in the global number of misses is found at K = 8 . For K > 8 , a little amount of useful information is added to the dictionary at the expense of increased coherence; (b) distance between recovered events (hits) and their corresponding simulated true event; (c) average amplitude of the events computed as hits, for noise levels σ = 0 (up), σ = 0.08 (middle) and σ = 0.12 (bottom). The bars indicate one standard deviation above and below the average. All simulated events have unity amplitude.
Figure 6. (a) percentage of misses (from 1000 simulated events) as a function of K, for three levels of noise, with OMPED running with a fixed number of five iterations (each of the 200 simulated acquisition had five scatterers). Each recovered scatterer distant more than 0.5 mm in any direction (axial or lateral) from the closest original simulated scatterer was computed as a miss. A minimum in the global number of misses is found at K = 8 . For K > 8 , a little amount of useful information is added to the dictionary at the expense of increased coherence; (b) distance between recovered events (hits) and their corresponding simulated true event; (c) average amplitude of the events computed as hits, for noise levels σ = 0 (up), σ = 0.08 (middle) and σ = 0.12 (bottom). The bars indicate one standard deviation above and below the average. All simulated events have unity amplitude.
Sensors 18 04097 g006
Figure 7. (a) histogram of final iteration (when e est e for the first time) for OMPED running with the SVD dictionary, for K varying from 2 to 10. Results from all values of K are summed. The total number of reconstructions is 5400. The 5th iteration was more frequently identified as final iteration, which is correct since all simulated acquisitions contained 5 scatterers; (b) example of evolution of e est and e along the iterations of OMPED. In this case, e est dropped below e at the 5th iteration, which was correctly identified as the final iteration. The simulated object contained five scatterers. White Gaussian noise with σ = 0.12 was added to the acquired data. OMPED was ran with the SVD dictionary with K = 8 .
Figure 7. (a) histogram of final iteration (when e est e for the first time) for OMPED running with the SVD dictionary, for K varying from 2 to 10. Results from all values of K are summed. The total number of reconstructions is 5400. The 5th iteration was more frequently identified as final iteration, which is correct since all simulated acquisitions contained 5 scatterers; (b) example of evolution of e est and e along the iterations of OMPED. In this case, e est dropped below e at the 5th iteration, which was correctly identified as the final iteration. The simulated object contained five scatterers. White Gaussian noise with σ = 0.12 was added to the acquired data. OMPED was ran with the SVD dictionary with K = 8 .
Sensors 18 04097 g007
Figure 8. Example of image simulated and reconstructed, from the dataset described in Section 6.1. The simulated data contains five scatterers and white Gaussian noise with σ = 0.08 . All images are normalized by the maximum absolute pixel value. (a) ground truth, with 5 unity amplitude scatterers randomly distributed over the ROI; (b) results from OMPED with the SVD dictionary ( K = 8 ). The algorithm correctly identified the 5th iteration as the final one; (c) results from OMP with regular model H . Seven iterations were run to show that, after the 4th iteration, the algorithm creates artifacts on the neighborhood of the leftmost scatterer instead of identifying the bottom-right scatterer present on the ground truth image. (d) Solution of the unregularized LS problem of Equation (4). The image is dominated by artifacts; (e) solution of the 1 -regularized problem of Equation (28). The penalization of the recovered amplitudes causes the suppression of most points on the resulting image. The chosen regularization parameter λ = 2.0691 minimizes the norm v v ^ , where v is the ground truth.
Figure 8. Example of image simulated and reconstructed, from the dataset described in Section 6.1. The simulated data contains five scatterers and white Gaussian noise with σ = 0.08 . All images are normalized by the maximum absolute pixel value. (a) ground truth, with 5 unity amplitude scatterers randomly distributed over the ROI; (b) results from OMPED with the SVD dictionary ( K = 8 ). The algorithm correctly identified the 5th iteration as the final one; (c) results from OMP with regular model H . Seven iterations were run to show that, after the 4th iteration, the algorithm creates artifacts on the neighborhood of the leftmost scatterer instead of identifying the bottom-right scatterer present on the ground truth image. (d) Solution of the unregularized LS problem of Equation (4). The image is dominated by artifacts; (e) solution of the 1 -regularized problem of Equation (28). The penalization of the recovered amplitudes causes the suppression of most points on the resulting image. The chosen regularization parameter λ = 2.0691 minimizes the norm v v ^ , where v is the ground truth.
Sensors 18 04097 g008

Share and Cite

MDPI and ACS Style

Rigo Passarin, T.A.; Wüst Zibetti, M.V.; Rodrigues Pipa, D. Sparse Ultrasound Imaging via Manifold Low-Rank Approximation and Non-Convex Greedy Pursuit. Sensors 2018, 18, 4097. https://doi.org/10.3390/s18124097

AMA Style

Rigo Passarin TA, Wüst Zibetti MV, Rodrigues Pipa D. Sparse Ultrasound Imaging via Manifold Low-Rank Approximation and Non-Convex Greedy Pursuit. Sensors. 2018; 18(12):4097. https://doi.org/10.3390/s18124097

Chicago/Turabian Style

Rigo Passarin, Thiago Alberto, Marcelo Victor Wüst Zibetti, and Daniel Rodrigues Pipa. 2018. "Sparse Ultrasound Imaging via Manifold Low-Rank Approximation and Non-Convex Greedy Pursuit" Sensors 18, no. 12: 4097. https://doi.org/10.3390/s18124097

APA Style

Rigo Passarin, T. A., Wüst Zibetti, M. V., & Rodrigues Pipa, D. (2018). Sparse Ultrasound Imaging via Manifold Low-Rank Approximation and Non-Convex Greedy Pursuit. Sensors, 18(12), 4097. https://doi.org/10.3390/s18124097

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop