CROSS-REFERENCE TO PRIOR APPLICATIONS

[0001] This patent application claims priority to U.S. Provisional Application Serial No. 62/425,336 filed on November 22, 2016. The application is hereby incorporated by reference herein in its entirety.

FIELD

[0002] This patent application pertains to ultrasound imaging of objects, especially in medical imaging, and more specifically relates to imaging at high spatial resolution that makes use of both magnitude and phase of echoes in an image reconstruction process that applies unique constraints to a fitting of echoes from the object of interest to echoes from an array of known scatter ers.

BACKGROUND

[0003] Super-resolution imaging or sub -wavelength imaging has become a popular domain of research particularly in medical imaging. Over the past couple of decades, the world of imaging has seen drastic improvement in the resolution as well as image quality with the help of improved hardware and efficient image processing algorithms. Ultrasound can be a preferred medical imaging modality due to its non-ionizing, non-invasive and relatively inexpensive nature. A desire to achieve high resolution in ultrasound imaging comparable to that of MRIs, CTs and X-ray tomography, has fueled extensive research aimed at improving resolution and speed of acquisition of ultrasound images. High frequency ultrasonic imaging systems often termed ultrasonic biomicroscope or UBM have become popular for imaging shallow tissue structures such as the skin and eyes as well as small animals. The UBM systems work at frequencies between 30-60 MHz and sometimes higher.

[0004] Over the past few years, compressive sensing [1] [2], deconvolution [3] [4] and adaptive/minimum-variance beamforming [5-9] have been researched and applied to the field of medical imaging. Compressive sensing, which is based on assuming sparsity at certain steps during the imaging process, has helped reduce the amount of data accrued and/or used for image reconstruction. It can be done in a few ways: a) By assuming a sparse scatterer distribution or b) assuming sparsity in the beamformed RF image or c) utilizing sparse RF data at the receiver [10-11]. The sparse scatterer map based compressive sensing methods discussed in [12] integrate single plane wave transmission for high frame rate imaging while achieving resolution comparable to that of conventional ultrasound. Another technique called Xampling [13-14], has imaged macroscopic perturbations from real cardiac ultrasound data with reduced speckle content. The authors in [15] and [ 16] discuss image reconstruction from sub-sampled raw RF data and present results that have comparable resolution to conventional ultrasound while using less than 30% of the original RF samples acquired before beamforming. Adaptive beamforming techniques have demonstrated improvement in image resolution with the help of reduced main lobe width and low side lobe levels. Adaptive beamformers calculate the weights for receive apodization based on the recorded data rather than using pre-calculated values. The authors in [5] and [6] report better resolution compared to delay-and- sum beamformers while using a smaller aperture and parallel receive beamforming by a minimum variance beamformer method. These and other techniques aim to achieve resolution equivalent to conventional ultrasound while reducing the RF data required for reconstruction. Super-resolution (SR) or submillimeter imaging, on the other hand, aims at reducing the dependence of resolution on pulse shape and width as well as the presence of speckles to achieve higher resolution.

[0005] There are super-resolution techniques that obtain images by combining many low resolution images and applying post image processing techniques. Super-resolution through time reversal acoustics can also be achieved due to the random nature of inhomogeneous media. Authors of [17] and [18] report results that indicate sub -wavelength imaging through time reversal acoustics. The phase-coherent multiple signal classification (MUSIC) method discussed in [19] is said to improve ultrasound resolution to a quarter of a wavelength. This method assumes a grid size smaller than the transducer length and the results are presented for frequencies in the 4-11 MHz range and at a depth of a couple centimeters. Several research groups [20-22] have used gas micro-bubbles for submillimeter ultrasound imaging in vascular systems that are in the near-field of imaging (few millimeters) and report tenfold improvement in resolution compared to conventional ultrasound. Another research group [23] reports improving the temporal resolution and acquisition times seen in ultrasound localized microscopy by applying a super-resolution optical fluctuation imaging method. Research presented in [24] discusses a blind 2-D deconvolution technique based on an improved phase-unwrapping technique applied to the pulse estimation and reports images with sharper tissue boundaries when compared to the images before deconvolution. Another group of authors [25] describes a methodology that consists of performing parametric modeling on the Fourier transform of the Hilbert transform of the RF data and achieves sub -wavelength resolution at higher frequencies (20MHz in their results). Clement et. al [26] describe a super-resolution image recovery technique based on Fourier spatial frequency spectrum analysis of the signals that are backprojected in the wave-vector domain to the focal plane. Their technique is used to estimate the size and location of the object and were able to detect a human hair with diameter 0 . 09 λ at 4.7 MHz at a few tens of millimeters depth.

[0006] Additional techniques for supper-resolution and related processing are discussed in [33]-[38], and some are understood to involve using modeling in which initial echoes from a model are processed with echoes from an actual object to ultimately generate an image of a region of interest in the object.

[0007] Most of the techniques referred to above aim for high resolution for superficial tissue imaging/characterization rather than deep tissue imaging at depths of a few hundred wavelengths.

SUMMARY

[0008] An object of the new approach described in this patent specification is super- resolution ultrasound imaging that involves the pre- and post-processing that estimates ultrasound images of an actual object from information derived from a priori knowledge of a sensing matrix. The desired images can be formed from the echoes (RF data) after only a single plane wave excitation of the object with the help of a least squares with novel constraints. The process can be extended to 3D.

[0009] One embodiment comprises an ultrasound transducer that sends ultrasound into an object from plural sending transducer elements and receives echoes from the object at fewer but no less than two elements of the transducer, without apodization or beam forming. A computer memory stores a sensing matrix of reflectance coefficients that when applied to an array of known scatterers produces echoes that said fewer receiving elements of the transducer receive or would have received from the scatterers in response to ultrasound that is the same as or approximates the ultrasound the transducer sends into the object. A computer processor estimates ultrasound properties of locations in the object that spatially relate to said array of scatterers by applying a bounded least squares estimation process to fit the echoes received from the object to a combination of the sensing matrix and properties of said location in the object. The estimation process constrains the properties of the locations in the object to positive values between zero and one. The computer then produces and displays an ultrasound image of the object as a function of said estimated ultrasound properties.

[0010] Another embodiment is an ultrasound imaging system comprising an ultrasound transducer having an array of transducer elements configured to detect both amplitude and phase of echoes from an object at fewer but not less than two of the transducer elements from which the transducer has sent ultrasound into be object to cause said echoes, a computer memory storing reflectance coefficients that when applied to an array of known scattering locations results in echoes that were or would have detected at said fewer transducer elements caused by ultrasound the transducer is configured to send into the object, and a computer processor configured to estimate ultrasound properties of locations in the object that spatially relate to said array of scatterers by applying a bounded least squares estimation process to fit the echoes received from the object to a combination of the sensing matrix and properties of said location in the object. The estimation process is configured to constrain the properties of said locations in the object to positive values between a minimum and a maximum. The computer processor is further configured to produce and display an ultrasound image of the object as a function of said estimated ultrasound properties. The constrained values can range from zero to one.

[0011] This patent specification further describes an ultrasound imaging method comprising sending ultrasound into an object from plural transducer elements and receiving echoes from the object at fewer but no less than two of the transducer elements, without apodization or beam forming, providing a sensing matrix that when applied to an array of known scatterers produces echoes that said fewer transducer elements receive or would have received from the scatterers in response to ultrasound approximating the ultrasound sent into the object, carrying out a bounded estimation process with a computer to fit the echoes received from the object to a combination of the sensing matrix and properties of locations in the object spatially related to said scatterers, wherein said estimation process is configured to constrain the properties of said locations in the object to positive values between zero and one, and producing and displaying an ultrasound image of the object as a computer-calculated function of said estimated ultrasound properties.

[0012] In another embodiment, an ultrasound imaging process comprises detecting both amplitude and phase of echoes from an object at fewer but not less than two of plural transducer elements from which ultrasound has been sent into the object to cause said echoes, providing a sensing matrix that when applied to an array of known scatterers produces echoes that said fewer transducer elements receive or would have received from the scatterers in response to ultrasound approximating the ultrasound sent into the object, carrying out an estimation process with a computer programmed to fit the echoes received from the object to a combination of the sensing matrix and properties of location in the object spatially related to said scatterers through a bounded least squares process, wherein said estimation process is configured to constrain the properties of said locations in the object to positive values between a minimum and a maximum, and producing and displaying an ultrasound image of the object as a computer-calculated function of said estimated ultrasound properties. The positive values can range from zero to one.

BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Fig. 1 illustrates an object, an ultrasound transducer, and an array of scaterrers superimposed on a region of interest on the object to provide a sensing matrix used in estimating an ultrasound image of a region of interest in the object.

[0014] Fig. 2 illustrates an example of creating a sensing or imaging matrix.

[0015] Fig. 3 illustrates an example of creating columns of sensing or imaging matrix A by compressing RF data matrix 2¾ by projecting the rows of Bij onto a random vector of positive l 's and negative l 's.

[0016] Figs. 4a-4c illustrate a comparison of the estimated images using two receive channels, where 4a shows the actual reflectivity pattern, 4b shows the estimated image using least squares without bounds, and 4c shows the estimated image using least squares with bounds.

[0017] Figs. 5a and 5 b show simulation results using a diagonal row of scatterers, where 5a is an actual reflectivity pattern and 5b is the estimated image using bounded least squares and two receiving channels. [0018] Figs. 6a and 6b are plots showing mean squared error vs (6a) frequency (6b) number of excitation cycles.

[0019] Fig. 7a shows an Axial Point Spread Function and Fig. 7b shows a Lateral Point Spread Finction.

[0020] Figs. 8a, 8b, and 8c illustrate experimental results for an object with two wires separated by 0.8 mm laterally, where 8a shows a conventional B-mode ultrasound image of the object with the wires on left and a magnification of the wire image at right, 8b and 8c show magnified images of the wires using the imaging technique described in this patent specification, at one location in 8b and at another location in 8c.

[0021] Figs. 9a, 9b, and 9c illustrate additional experimental results for an object with two wires separated by 0.8 mm laterally and axially, where 9a shows a conventional ultrasound image of the object with the wires on left and a magnification of the wire image at right, 9b shows an interpolated image of the two wires separated by 0.8 mm using the new approach described in this patent specification, and 9c shows the wires at a different location

[0022] Figs. 10a, 10b, 10c, and lOd shows results of an experiment where fishing wires are placed randomly: 10a shows a magnification of a conventional ultrasound image with three wires, 10b shows an image of the three wires produced as described in this patent specification, 10c shows a magnification of four wires in a conventional ultrasound image, and lOd shows an image of the four wires produced with the new approach of this patent specification.

[0023] Figs. 11a and 1 lb illustrate further experimental results: 1 la shows a conventional ultrasound image of two wires separated by 0.8 mm passing through beefsteak (full image on left and magnification on right) and Fig. 1 lb shows an image of the same object produced with the new approach described in this patent specification.

[0024] Fig. 12a shows further experimental results using the new approach described in this patent specification on the same object (two wires separated by 0.8 mm in a beefsteak): 12a shows a single frame image and Fig. 12b shows an averaged image, over several frames.

[0025] Fig. 13 illustrates main steps of an example of estimating ultrasound properties of locations in an object.

[0026] Fig 14 is a block diagram of a system for estimating ultrasound properties of locations in an object. DETAILED DESCRIPTION

[0027] A detailed description of examples of preferred embodiments is provided below. While several embodiments are described, the new subject matter described in this patent specification is not limited to any one embodiment or combination of embodiments described herein, but instead encompasses numerous alternatives, modifications, and equivalents. In addition, while numerous specific details are set forth in the following description to provide a thorough understanding, some embodiments can be practiced without some of these details and even without all of the described details. Moreover, for clarity and conciseness, certain technical material that is known in the related technology have not been fully described in detail, to avoid unnecessarily obscuring the new subject matter described herein. It should be clear that individual features of one or several of the specific embodiments described herein can be used in combination with features of other described embodiments or with other features. Further, like reference numbers and designations in the various drawings indicate like elements.

[0028] Referring to Fig. 1, an ultrasound transducer 102 has plural transducer elements such as 48 elements and is acoustically coupled with an object or body 104 and configured to send transmitted ultrasound such as a plane wave into object 104 from all the transducer elements, without apodization of beam forming. Transducer 102 detects echoes from object 104 at fewer but not less than 2 of its elements, for example at the two end elements in the case of a linear array of 48 elements. The transducer can have a different, typically greater, number of transducer elements, and can detect echoes at more than two of the elements so long as the detecting elements are fewer, and typically much fewer, than the elements that transmit ultrasound into object 104. And, the transducer elements can be arranged other than in a linear array, for example in plural linear arrays, in one or more curved arrays, or in a 2D pattern that can be periodic or non-periodic. Also illustrated in Fig. 1 is a grid 106 defining an array of scatterers or scattering points 106a spaced from each other in one direction by distances dx and by distances dz in another direction. Distances dx and dz can be the same or different. Grid 106 can be square as illustrated, or rectangular, or in some other shape, and can be periodic as illustrated or non- periodic, and can be defined in a coordinate system other than an x-z system, such as in polar coordinates. Grid 106 is a notional rather than an actual grid and, as described below, is used in estimating an ultrasound image representing locations in object 104 that spatially match locations of scatterers 106a when grid 106 is notionally superimposed on a region of interest in object 104. [0029] In principle, the new approach described in this patent specification derives or provides a sensing matrix that, when applied to a known array of scatterers, results in echoes that those scatterers produce or would have produced in response to ultrasound that approximates the ultrasound that transducer 102 sends into object 104. A process with unique bounding fits the echoes received from the object to a combination such a product of the sensing matrix and ultrasound properties of locations in the object that spatially relate to said scatterers, and produces an ultrasound image of a region of interest that comprises said locations in the object. The process can include bounded least squares operation that constrains the properties being estimated to a range of positive values such as zero to one. An image of a region of interest in the object is produced as a function of the estimated ultrasound properties.

[0030] First, a theoretical basis is laid out below, followed by implementation examples.

[0031] A few essential equations are described below. A fluid model is assumed for propagation of ultrasound waves, so the waves that travel in a medium used in this discussion are longitudinal in nature and have an associated wavenumber k = ω/c, where ω = 2nf and /

du

is the spatial frequency. Particle velocity in the medium is given by v =— , where u is the particle displacement. For convenience, the particle velocity is also expressed as the gradient of velocity potential, V0. The pressure is then given by p =— p ^.

[0032] The plane wave equation (three dimensional) that governs acoustic propagation in an ideal medium is given by [27]:

1 d ^z

4>(x, y, z, t) -—— (x, y, z, t) = 0 (1) where c is the propagation velocity in the medium, φ (x, y, z, t) is the pressure or the velocity potential at the location (x, y, z). The Helmholtz wave equation in (1) is in general assumed to be linear and time shift invariant in both time and space. In the frequency domain, Eq. (1) can be expressed as, 2φ + /c ² <p = 0 (2) [0033] where φ is the Fourier transform of <p. k is the wavevector and can be broken down into its projections in the x, y, and z directions as k ² = k _% + k _y + k . Analogous to the time domain waveform having a spectrum consisting of a collection of frequencies, the acoustic field of a transducer has an angular spectrum of plane waves [28]. As such, a Fourier transform relation can be established between the amplitude of the source and the spatial frequency distribution as discussed in Szabo [29].

[0034] The finer details of an object are often recovered from the higher spatial frequencies. In a way, the imaging system acts as a filter where the resulting signal is a convolution of the object, g, with the point spread function of acoustic system, say h. A mathematical model describing the imaging process is shown in Eq. (3 ). y = h(x) * g(x) (3)

In the frequency domain, this equation becomes,

Y = H(f)G(f) (4) where H is called the modulation transfer function (MTF), G is the Fourier transform of g and

Y is Fourier transform of_y. Due to the finite aperture size of the transducer, there is a low-pass cutoff frequency above which information is typically lost. Typically, the MTF has spectral notches (or zeros) rendering the convolution operation in Eq. (3) non-invertible. Consequently, only the spectral components below the first spectral notch (i.e., spatial frequency components below the diffraction limit) can be uniquely recovered. During the measurement process in most imaging systems, only coarse features of the object g are obtained, thus not being able to resolve the finer details. Per results presented in [30] by Hunt, when a wave is reflected by a scatterer, it introduces high frequency spatial content which affects the spatial frequencies below the cutoff spatial frequency o/c. As such the problem of super-resolution is one that involves recovering the fine details of an object being imaged with the coarse measurements from the spectrum below the cut-off spatial frequency. [0035] Work reported in [30-31] informs that super-resolution can be achieved if certain conditions like positivity, compactness, or sparsity are assumed for the support of the signal g during the process of image reconstruction. In some cases, this has been achieved with the help of a priori information [12], [26], and [31]. This can be recognized as a starting point for the sub-wavelength ultrasound imaging technique described in this patent specification.

[0036] A discussion follows of main steps of an example of image reconstruction per principles described in this patent specification.

[0037]

[0038] The analysis of ultrasound systems [29] begins with the pulse-echo equation of the echo, r(t), as given in Eq. (5). ^r ^ ⁼ ^^ -C _o -C _o ^^ _' ^ _' ² ) ⁷¹ ^ - 2c- ¹ z)q ² (x, y, z)dxdydz (5) where q ² (x, y, z) = zq x, y, z). q x, y, z) is the transducer field pattern, c is th e speed of sound in the medium, and k is the wavenumber. μ _α accounts for the attenuation in the medium being imaged. R(x, y, z) is the scatterer strength or reflectivity. In traditional B-mode ultrasound, only the magnitude of the received signal R(x, y, z) is used during image reconstruction to form each A-line. But in the super resolution ultrasound system described in this patent specification, both the phase and amplitude are used for the image reconstruction process. In the technique example described below, a plane wave is sent out from a linear transducer array and the received echoes are stored without applying apodization or beamforming and thus being able to conserve the phase information.

[0039] Consider the case of two-dimensional imaging in the x-z plane where the axis of propagation is in the z direction (refer to Fig. 1). A tran sducer 1 02 send s a pl ane wave ultras ound i nto an object 104 (a medium such as tissue). A two-dimensional grid 106 o f s c att e r e r s 1 06 a is formulated to overlay the object with Nx and Nz as the dimensions in the x and z directions, respectively. We then define a reflectivity matrix #« _xX « _z , where each element y corresponds to a reflectivity of a point scatterer 106a in the object at every grid point x, and zj. A column vector ¾ _xWzXl is then formed by concatenating the columns of the matrix R ⁹ as shown below, where R^is the j ^th column of the matrix. XN _x N _z xl

[0040] Now we form a sensing or imaging matrix A. The matrix A can be described as a set of column vectors that correspond to the spatial impulse response (or signal received by the transducer array, which can be called echoes detected or received by transducer 102) when only one point scatterer 106a with maximum reflectivity is present in the object as shown in Eq. (6). Let us assume that the array has Nc receive channels and each channel takes a total of Ns samples per frame per acquisition. The imaging matrix A is then of the size Nc A¾ by Nx Nz. For example, a column vector Ajj _cNsXl corresponds to the signal samples (echoes) received by the transducer array 102 from a point scatterer 106a located at (1,1) on the grid of size (Nx, Nz).

[0041] Figure 2 illustrates the formation of one column vector of sensing or imaging matrix A. Since matrix A is calculated prior to the process of imaging an object 104, it can be considered a priori information. The sensing or imaging matrix can be written as below,

A = \A ¹ A ² . .. A ^N * ^N A (6)

The received RF data (echoes) from grid points of any target phantom will be a linear combination of the vectors in the imaging matrix A. From here the imaging model can be formulated as shown below,

¾ _c W _s xl ⁼ ^N _c N _s xN _x N _z ^N _x N _z xl + EN _C N _S X.1 )

[0042] where E is a column vector of samples of white Gaussian noise process with covariance matrix CE, F is the measured raw RF data, and X is a column vector of reflectance coefficients to be estimated. Estimation in the ultrasound imaging method in this example is carried out using the least squares estimation algorithm, which finds an estimate of the column vector X _{N N xl} that minimizes L2-norm of the error term ε = AX— Y as shown in Eq. (8).

¾ _xWz xi = argmin(||^ - r|||) (8) x Even though the unbounded least-squares estimation method in Eq. (8) may be computationally efficient as compared to other methods (e.g., sparsity enhancing estimation methods), it may prove non-robust under certain conditions when the imaged object deviates from the assumed imaging model in Eq. (7). Examples of such deviations are scatterers outside the region of interest (or scatterers outside the grid), point scatterers not aligned with grid intersections, reverberations and dispersion in the imaged object, and non-stationary noise. To improve the performance of the estimation algorithm in (8), a weighted least-squares may be used instead as shown in Eq. (9), where the less accurate observations ^/ of the received vector F are weighted to produce smaller effect on the estimation of X.

XN _X N _Z XI = argmin(∑" " ^s w _fc £¾ , 0 < w _k ≤ 1 (9) x

[0043] Another way to improve the performance of the estimation method in (8) is to use a bounded least-squares estimation method where the estimates of the reflectance coefficients are bounded to some region [a,b] as shown in Eq. (10). For example, the estimated reflectance coefficients may be bounded to positive values from zero to one.

¾ _xWz xi = argmin(||^ - r|||) (10)

Xe[a,b]

[0044] The discussion below and Fig. 4 illustrate that bounded least squares method is a preferred technique for estimation of the image coefficients in X when compared to the case of unbounded least squares.

[0045] Lastly, further improvement in the reflectance coefficients estimation may be achieved by using an additional LI -norm in the cost function as shown in Eq. (1 1), where the weighting coefficient λ is a regulation parameter. This estimation method takes advantage of sparsity in scatterer distribution in typical ultrasound images and is known as LASSO estimation method, [39]. argminai^ - ni + imii) (11)

[0046] The sensing or imaging matrix A is typically bandlimited (or sparse in 2-D frequency domain) and can be 'compressed' to a smaller rank matrix without losing too much of the information available in RF data (echoes). Therefore, in another embodiment described in this patent specification, the sensing or imaging matrix A can be obtained by compressing the RF data (echoes) from Nc transducer channels each taking Ns samples to a fewer number of samples as described below. For each point in the grid at the location (x, , zj) the RF signal from Nc transducer elements of the ultrasound probe is either pre-calculated or experimentally determined forming a matrix 2¾ of the size Nc by Ns, where Ns is the total number of samples received by one transducer element (note that each grid point will have different matrix B unique to that grid point). The corresponding column i4 ^fc of the imaging matrix A is then formed by compressing matrix Bij with the use of a linear operator L as shown in Eq. (12), where the column vector A ^k is of size N 1 ( N≤N _C N _S ).

A xi = £{BN _C XN _S } (1 )

[0047] This reduces the size of the estimation problem and associated computational costs. An example of the RF data compression and formation of one column vector of the compressed imaging matrix A is shown in Fig. 3. As shown in the figure, the rows of the RF data matrix Bi _j are projected onto a pseudo-random vector of I 's and -I 's of size N _c x 1 (i.e., each sample akm of the column vector i4 ^fc is calculated as an inner product between an m" ¹ row of Bi _j and m" ¹ pseudo-random vector of I 's and - I 's). Note that the pseudo-random vectors can be different for each row of Bij

[0048] Simulation results described below have confirmed the theoretical explanation given above for the new approach to super-resolution ultrasound imaging. Simulation were set up with the help of the open source Acoustics toolbox from k-wave [32] along with MATLAB

(Mathworks Inc, Natick MA). A two-dimensional (x-z) grid was set up in a simulated object or medium with properties similar to homogeneous tissue (c = 1540 m/s and density = 1000 kg/m ³ ). The simulated scatterers were placed at a depth close to 5 cm. The simul ated excitation signal consisted of a 5-cycle sinusoid with center frequency 1.875 MHz, windowed by a Gaussian profile to mimic the signal emitted from a transducer such as 102 (Fig. 1) upon excitation. In this example the technique uses plane wave excitation, and no apodization or focusing was applied on either transmit or receive. The acoustic source (such as transducer 102 in Fig. 1) used for simulations was assumed to have 48 active elements on transmit and 2 active elements on receive. Using one channel alone may be sufficient for imaging in 1-D but for 2-D imaging, at least two receive channels (such as for example the two end elements of transducer 102) are desirable so as to provide accurate estimations in the presence of system noise. In case of 2-D imaging, the signal received by one element of the array of transducer elements from a scatterer such as 106a at one depth may appear the same as the signal received from a different scatterer at another point that lies on the arc intersecting the first point. As such, a better estimate will be provided if at least two transducer channels (elements of transducer 102) are used to track the point scatterer 106a in the medium.

In the simulation examples, firstly the impulse response of each point in a grid of size Nx = 11 and Nz = 5 was stored in computer memory. For simplicity, uniform grid spacing of 0.4 mm was considered even though the algorithm used in the simulation does not require uniform grid spacing in both dimensions. Next, an object (phantom) consisting of scatterers of higher density and acoustic speed than the surrounding medium was excited using the aforementioned transducer and excitation signal. The phantoms were constructed as disc structures with an acoustic speed of 3000 m/s. After receiving the echoes (RF data), image reconstruction was performed using the bounded least squares estimation equation shown in Eq. (10) with bounds from zero to one.

[0049] The importance of using least squares with bounds is shown in Fig. 4, where the figures show a comparison between unbounded and bounded least squares estimation. The columns of the imaging matrix A are highly correlated and as a result any response from the phantom that is not perfectly aligned with one of the columns of the matrix A will have nonzero projections onto more than one column vector of the matrix A. For example, if the scattering point of the phantom is located off the grid, its response will project not only onto the impulse response from the nearest neighbors, but rather onto impulse responses from a large neighborhood. As a result, similar to the linear deconvolution methods the decoded image with unbounded least-squares estimation exhibits widely spread point spread function indicating poor resolution.

[0050] Super-resolution can be achieved if the conditions of positivity and compactness are met [30], which is why the images estimated using bounded least squares achieve super- resolution. Additional simulations were run for increasing number of point scatterers 106a and various patterns confirming that least squares estimation with bounds as described in this patent specification provides a good estimate of the target phantom (object) provided at least two channels (elements of transducer 102) are used. Figs. 5a and 5b show the estimated image of a diagonal line of scatterers (Fig. 5(b)) in comparison to the actual reflectivity pattern (Fig. 5(a)). Here two receive channels were used for decoding.

[0051] Further simulations were run with an increasing number of receiving channels. In the absence of noise, adding more than two channels did not provide substantial improvement to image reconstruction quality for the setup used in the simulations. Still other simulations were run to investigate the effects of grid size and change in frequency on the image

reconstruction. The results shown in Fig.6(a) indicate that increasing the frequency for a grid size of 0.4 mm at a given depth of 4.86 cm decreases the mean squared error (MSE) of the estimated images. At frequency of 8 MHz and 16 MHz, the error starts to flatten out since attenuation starts becoming prominent. From simulations, we see that excitation frequency and grid spacing (dx, dz) follow the same error trend. For example, reducing the frequency by a factor of two is similar to reducing the grid size by a factor of two. As shown in Fig. 6(b), the MSE also increases with increasing number of cycles in the excitation signal which goes hand in hand with the fact that resolution decreases as the spatial pulse length increases where spatial pulse length is the product of number of cycles and excitation wavelength.

[0052] Figures 7a and 7b shows the point spread function(PSF) of the simulated ultrasound system by sweeping over a point in the axial direction. The axial PSF has a full width half maximum (FWHM) of 0.1264 mm at the excitation frequency of 1.865 MHz. The lateral PSF has a wider profile as expected, indicating lower resolution compared to axial resolution. The lateral PSF has a FWHM of 0.436 mm.

[0053] Experimental results further confirmed the new approach described in this patent specification. Experiments were performed using Verasonics VI ultrasound scanner (Verasonics, Inc. Kirkland WA 98034) connected to a 96-channel phased array ATL probe, P4-1 (ATL Ultrasound, Inc., Bothel WA 98041). Only the first half of the aperture (48 channels) was used for transmission while two channels from the same aperture are used on receive. Only the first half of the transducer aperture (14.16 mm of total aperture size of 28 mm) was used on transmission, to ensure imaging in the far field region while maintaining a sufficiently high signal to noise ratio (S R). As in the simulations, only two channels were used on receive to uniquely recover the 2-D target. The ultrasound phantom consisted of fishing wires in a tank filled with degassed water. The fishing wires are made of nylon and are of 0.2 mm in diameter. The probe was placed perpendicular to and roughly 8.9-9.2 cm above the fishing wires. The excitation frequency was set to 1.875 MHz corresponding to a wavelength of 0.833 mm for acoustic speed of 1540 m/s. Raw RF data (echoes) was collected and then passed through a FIR bandpass filter in MATLAB which removes noise outside the bandwidth of interest. This was followed by the image reconstruction algorithm. No apodization, focusing and beamforming were applied. The grid size used for the experiments was similar to the dimensions in the simulation setup and was large enough to cover the phantoms used for the experiments.

[0054] The experiments used were done in two-dimensions i.e. the x-z plane as illustrated in Fig. 1, as extension to 3 dimensions can be achieved by extending grid 106 (in Fig. l) to 3D space. The sensing or imaging matrix A was populated using a simple Nz=5, Nx=l l grid spanning 2 mm axially and 4.4 mm laterally, where each grid intersection represented a point scatterer. The grid points were taken 0.4 mm apart in the axial and lateral directions.

[0055] The first two-dimensional image involved two fishing wires separated by 0.8 mm laterally. The results are shown in Figs. 8a, 8b, and 8c, where 8a shows a conventional B-mode ultrasound image of the object with the wires on left and a magnification of the wire image at right, and 8b and 8c show magnified images of the wires using the imaging technique described in this patent specification, at one location in 8b and at another location in 8c. Some background noise is visible in the images and can be reduced by averaging the reconstructed images. Note that the axes of the magnifications of the conventional ultrasound images shown in the figures are units of wavelength.

[0056] Figs. 9a, 9b, and 9c illustrate additional experimental results for an object with two wires separated by 0.8 mm laterally and axially, where 9a shows a conventional ultrasound image of the object with the wires on left and a magnification of the wire image at right, 9b shows an interpolated image of the two wires separated by 0.8 mm using the new approach described in this patent specification, and 9c shows the wires at a different location.

[0057] Figs. 10a, 10b, 10c, and lOd shows results of an experiment where fishing wires are placed randomly: 10a shows a magnification of a conventional ultrasound image with three wires, 10b shows an image of the three wires produced as described in this patent specification, 10c shows a magnification of four wires in a conventional ultrasound image, and lOd shows an image of the four wires produced with the new approach of this patent specification. Additional imaging was done of three and four wire phantoms using a 9 x 11 (Nz = 9, Nx = 11) grid which spanned 3.6 mm axially and 4.4 mm laterally. The area outlined in the box in Fig.10a and 10c approximates the area shown in Fig.10 (b) and (d) respectively.

[0058] Figs. 11a and 1 lb illustrate further experimental results: 1 la shows a conventional ultrasound image of two wires separated by 0.8 mm passing through beefsteak (full image on left and magnification on right) and Fig. 1 lb shows an image of the same object produced with the new approach described in this patent specification. The speed of sound in water is roughly 1540 m/s and has an acoustic impedance of 1.483 g/cm ² sec x lO ⁵ and density equal to 1 g/cm ³ . Acoustic waves get reflected or refracted at boundaries between objects with different acoustic properties. For a real-life ultrasound imaging scenario, an experiment was performed to see the resolution capability of the technique behind living tissue. This was performed using a piece of steak of 1 cm thickness. Beef has an acoustic impedance of 1 68 g/cm ² sec x lO ⁵ and density of 1.08 g/cm ³ . As can be seen in Fig. 1 lb, the signal is attenuated due to the tissue layer above the wires while the resolution is clearly maintained.

[0059] Fig. 12a shows further experimental results using the new approach described in this patent specification on the same object (two wires separated by 0.8 mm behind a beefsteak): 12a shows a single frame image and Fig. 12b shows an averaged image, over several frames. The images presented before Figs. 12a and 12b were not averaged. In order to reduce the noise present in the estimated images, averaging over multiple frames has proven to help. Figure 12b shows the estimated image from the experiment with steak after 15 averages in comparison to a single frame of the image.

[0060] Fig. 13 illustrates main step of a process for estimating ultrasound properties of locations in an object according to the above description. In step 130 the process selects an array of known scatterers. In a simplified example, the known scatterers can be high-density point scatterers in water or fat. In step 132, the process generates a sensing matrix by calculations, with or without deriving actual echoes of the known scatterers by ultrasound measurements. In stepl34, the process sends ultrasound into the object and receives echoes as described above, and in step 136 fits the received echoes to a combination such as a product of the sensing matrix and a vector of ultrasound properties of locations in the object. The fitting process in step 136 may be performed by finding ultrasound properties of locations in the object that minimize fitting error in a least-square sense.

[0061] Fig. 14 illustrates a block diagram form of a system for carrying out the estimation of ultrasound properties of locations in an object using the processes described above. Ultrasound transducer 102 sends ultrasound into object 104 and receives echoes from the object as described above. An ultrasound engine 1402, which can be the engine of a conventional processor such as in the examples given in the background section of this patent specification, controls transducer 102 and processes the echoes received from the transducer by carrying out the processes described above, using programming that a person of ordinary skill in programming can implement based on the description and equations explained in this patent specification.

Computer storage 1404 stores the sensing matrix (based on the known array of scatterers).

Computer storage 1406 stores the initial values of ultrasound properties of locations in the object before the fitting process and their modified values that result from fitting echoes from the object, which are stored in computer storage 1408, to the combination such as a product of the sensing matrix and a vector of ultrasound properties of locations in the object. From the above- described fitting, ultrasound engine 1402 estimates ultrasound properties of locations in object 104 and displays them as an image or in another form on image display 1410.

[0062] Although the foregoing has been described in some detail for purposes of clarity, it will be apparent that certain changes and modifications may be made without departing from the principles thereof. It should be noted that there are many alternative ways of implementing both the processes and apparatuses described herein. Accordingly, the present embodiments are to be considered as illustrative and not restrictive, and the body of work described herein is not to be limited to the details given herein, which may be modified within the scope and equivalents of the appended claims. References:

[1]. E. J. Cande and M. B. Wakin, "An introduction to compressive sampling," Signal

Processing Magazine, IEEE, vol. 25, no. 2, pp. 21-30, 2008.

[2]. E. J. Candes and T. Tao, "Decoding by linear programming," Information Theory, IEEE

Transactions on, vol. 51, no. 12, pp. 4203-4215, 2005.

[3]. O. Michailovich and D. Adam, "Phase unwrapping for 2-d blind deconvolution of

ultrasound images," Medical Imaging, IEEE Transactions on, vol. 23, no. 1, pp. 7-25,

2004.

[4]. T. Taxt and G. V. Frolova, "Noise robust one-dimensional blind deconvolution of

medical ultrasound images," Ultrasonics, Ferroelectrics, and Frequency Control, IEEE

Transactions on, vol. 46, no. 2, pp. 291-299, 1999.

[5]. J.-F. Synnevag, A. Austeng, and S. Holm, "Benefits of minimum -variance beamforming in medical ultrasound imaging," Ultrasonics, Ferroelectrics, and Frequency Control,

IEEE Transactions on, vol. 56, no. 9, pp. 1868-1879, 2009.

[6]. S. Holm, J. Synnevag, and A. Austeng, "Capon beamforming for active ultrasound

imaging systems," in Proc. IEEE, 13th DSP Workshop, 2009.

[7]. J. A. Mann and W. Walker, "A constrained adaptive beamformer for medical ultrasound:

Initial results," in Ultrasonics Symposium, 2002. Proceedings. 2002 IEEE, vol. 2.

IEEE, 2002, pp. 1807-1810.

[8]. I. K. Holfort, F. Gran, and J. A. Jensen, "P2b-12 minimum variance beamforming for high frame-rate ultrasound imaging," in Ultrasonics Symposium, 2007. IEEE. IEEE,

2007, pp. 1541-1544.

[9]. B. M. Asl and A. Mahloojifar, "Eigenspace-based minimum variance beamforming

applied to medical ultrasound imaging," Ultrasonics, Ferroelectrics, and Frequency Control, IEEE Transactions on, vol. 57, no. 11, pp. 2381-2390, 2010.

[10]. C. Quinsac, A. Basarab, J. Gregoire, and D. Kouame, "3d compressed sensing ultrasound imaging," in Proc. IEEE Int. Ultrason. Symp.(IUS), San Diego, USA, 2010. [11]. C. Quinsac, A. Basarab, J.-M. Girault, and D. Kouame, "Compressed sensing of ultrasound images: sampling of spatial and frequency domains," in Signal Processing

Systems (SIPS), 2010 IEEE Workshop on. IEEE, 2010, pp. 231-236.

[12]. M. Schiffner, T. Jansen, and G. Schmitz, "Compressed sensing for fast image acquisition in pulse-echo ultrasound," Biomedical Engineer- ing/Biomedizinische

Technik, vol. 57, no. SI-1 Track-B, pp. 192-195, 2012.

[13]. N. Wagner, Y. C. Eldar, A. Feuer, G. Danin, and Z. Friedman, "Xampling in ultrasound imaging," CoRR, vol. abs/1104.5327, 2011.

[14]. N. Wagner, Y. C. Eldar, A. Feuer, and Z. Friedman, "Compressed beamforming applied to b-mode ultrasound imaging," in Biomedical Imaging (ISBI), 2012 9th IEEE

International Symposium on. IEEE, 2012, pp. 1080-1083.

[15]. D. Friboulet, H. Liebgott, and R. Prost, "Compressive sensing for raw rf signals

reconstruction in ultrasound," in Ultrasonics Symposium (IUS), 2010 IEEE. IEEE, 2010, pp. 367-370.

[16]. H. Liebgott, R. Prost, and D. Friboulet, "Pre-beamformed rf signal reconstruction in

medical ultrasound using compressive sensing," Ultrasonics, vol. 53, no. 2, pp. 525-533, 2013.

[17]. P. Blomgren, G. Papanicolaou, and H. Zhao, "Super-resolution in time-reversal

acoustics," The Journal of the Acoustical Society of America, vol. I l l, no. 1, pp. 230- 248, 2002.

[18]. A. J. Devaney, " 1 super-resolution processing of multi-static data using time reversal and music," 2000.

[19]. Y. Labyed and L. Huang, "Super-resolution ultrasound imaging using a phase-coherent music method with compensation for the phase response of transducer elements," Ultrasonics, Ferroelectrics, and Frequency Control, IEEE Transactions on, vol. 60, no. 6, pp. 1048-1060, 2013.

[20]. M. A. O'Reilly and K. Hynynen, "A super-resolution ultrasound method for brain

vascular mapping," Medical physics, vol. 40, no. 11, p. 110701, 2013.

[21]. B. Cox and P. Beard, "Imaging techniques: Super-resolution ultrasound," Nature, vol.

527, no. 7579, pp. 451-452, 2015. [22]. C. Errico, J. Pierre, S. Pezet, Y. Desailly, Z. Lenkei, O. Couture, and M. Tanter,

"Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging," Nature, vol. 527, no. 7579, pp. 499-502, 2015.

[23]. T. Dertinger, R. Colyer, R. Vogel, J. Enderlein, and S. Weiss, "Achieving increased resolution and more pixels with superresolution optical fluctuation imaging (sofi),"

Optics express, vol. 18, no. 18, pp. 18 875-18 885, 2010.

[24]. T. Taxt and R. Jirik, "Superresolution of ultrasound images using the first and second harmonic signal," Ultrasonics, Ferroelectrics, and Frequency Control, IEEE

Transactions on, vol. 51, no. 2, pp. 163-175, 2004.

[25]. D. Kouame and M. Ploquin, "Super-resolution in medical imaging: An illustrative

approach through ultrasound," in Biomedical Imaging: From Nano to Macro, 2009.

ISBF09. IEEE International Symposium on. IEEE, 2009, pp. 249-252.

[26]. G. Clement, J. Huttunen, and K. Hynynen, "Superresolution ultrasound imaging using back-projected reconstruction," The Journal of the Acoustical Society of America, vol.

118, no. 6, pp. 3953-3960, 2005.

[27]. J. N. Wright, "Image formation in diagnostic ultrasound," 1997.

[28]. J. W. Goodman, "Introduction to fourier optics," 2005.

[29]. T. Szabo, Diagnostic Ultrasound Imaging: Inside Out, ser. Academic Press series in biomedical engineering. Elsevier Academic Press, 2004. [Online]. Available:

[30]. B. R. Hunt, "Super-resolution of imagery: understanding the basis for recovery of spatial frequencies beyond the diffraction limit," in Information, Decision and Control, 1999. IDC 99. Proceedings. 1999. IEEE, 1999, pp. 243-248.

[31]. E. J. Candes and C. Fernandez-Granda, "Towards a mathematical theory of super- resolution," CoRR, vol. abs/1203.5871, 2012.

[32]. B. E. Treeby, J. Jaros, A. P. Rendell, and B. Cox, "Modeling nonlinear ultrasound

propagation in heterogeneous media with power law absorption using a k-space pseudospectral method," The Journal of the Acoustical Society of America, vol. 131, no. 6, pp. 4324-4336, 2012. [33]. F.S. Foster, C.J. Pavlin, K.A. Harasieqicz, D.A. Christopher and D.H. Turnbul,

"Advances in Ultrasound Biomicroscopy," Ultrasound in Med. & Biol., Vol. 26, No. 1, pp. 1-27, 2000.

[34]. F. Viola, M.A. Ellis and W.F. Walker, "Time-Domain Optimized Estimator for

Ultrasound Imaging: Initial Development and Results," IEEE Trans Med Imaging. 2008

January ; 27(1); 99-110. Doi 10.1109/TMI. 2007.93579.

[35]. M.A. Ellis and W.F. Walker, "Super-Resolution Image Reconstruction with Reduced

Computational Complexity," 2009 IEEE International Ultrasound Symposium

Proceeding, pp. 2351-2354.

[36]. M.A. Ellis and W.F. Walker, "Super-Resolution Image Reconstruction with Diffuse

Source Models," Ultrasound Med Biol. 2010 June ; 36(6): 967-977.

Doi: 10.1016/j .ultrasmedbio.2010.03.002..

[37]. Y. Yankelevsky, Z. Freedman and AFeuer, "Component Based Modeling of Ultrasound

Signals," arXiv: 1603.00273vl .

[38]. W.F. Walker and M. Ellis, U.S. Patent No. 8,818,064 B2, August 26, 2014.

[39]. R. Tibshirani, "Regression shrinkage and selection via the lasso," J. Roy. Stat. Soc. B, vol. 58, pp. 267-288, 1996.