Rician Beamforming: Despeckle Method via Coarray Projection Stochastic Analysis

Current computation capabilities normal, Delay and Sum (DAS) and Total Focusing Method (TFM) allow new definitions for beamformers. In this paper, a new beamforming method is proposed. It exploits diversity across pixel data after focusing process. The method is based on statistical analysis and sparse array concept. It avoids common aberrations introduced by beamforming process without loosing the original image texture, producing a better quality image. We evaluate the proposed method through simulation on standard phantoms. Images resulting from our novel method display better quality and provide localised estimations on image noise.


Introduction
Ultrasonic imaging offers numerous advantages over other types of medical imaging, specifically it is low cost, harmless and an ideal alternative for outpatient basis. However, compared to higher-cost or non-innocuous medical imaging techniques, ultrasonic imaging presents lower contrast ratio compromising interpretation of output images. This is because of the physical characteristics of tissues under study, due to low acoustic impedance differences [1]. Other considerations relate to the nature of signals used in this medical imaging modality and ultrasonic interactions, developing into the so-called speckle phenomenon. This introduces a granular aberration, due to constructive and destructive sums within tissue structures. Leading to contrast ratio lowering and aberrations depending on the point of view of the aperture and the structure geometry itself [2]. These hinder the processing of ultrasonic images from segmentation [3] and registering [4], required for Computer Aided Detection (CAD) processing chains [5][6][7].
In recent years, different strategies to lower the influence of this phenomenon have been published. Basically, these strategies can be grouped into three categories: based on spatial compounding [8], based on frequency compounding [9,10] and those based on post-filtering processes on image domain [11][12][13]. None of these methods completely mitigate the aberrations present in the image, but produce notably better images.
In the proposed method two different noise types are considered: Gaussian additive noise produced by heat and electronic noise due to transducers and the instrumentation, as well as backscattering noise. The latter could be modelled as a multiplicative noise, though produced aberrations are found to be convolutional [11].
In this contribution, we report on a new approach based on Full Matrix Capture (FMC), TFM techniques-for capturing and focusing-and sparse random arrays processing. Sensibility improvement to transduction systems [14] and ultrasound frontend massive data management have contributed to the popularity of these chained techniques in recent years. Image quality can be improved while increasing complexity of acquisition systems. This is not only due to this cascade of techniques, but to the ability to extract disaggregated data related to each transmitter-receiver pair. Additionally, after applying the Synthetic Aperture Focusing Technique (SAFT), disaggregated data of each point within the region of interest (ROI) can also be extracted. This allows for new strategies to mitigate the effect of secondary lobes and low directivity of the main lobe that translate into speckle noise in the resulting ultrasonic images [15].
In light of this, we describe a new postprocessing approach: we propose to use virtual sparse arrays to exploit spatial diversity on disaggregated data from SAFT. In this way, robust estimates can be obtained, beyond standard limitations of traditional beamformers lower resolution and speckle effect [8]. The resulting stochastic beamformer produces better contrasted images with better defined tissue interfaces, as well as better estimates on structures spatial reflectivity.
This new methodology contributes to a reliable noise characterisation at specific positions in the analysis area. Modern speckle reduction techniques based on homomorphic filtering [16], diffusion models [17], or total variation estimation [18], could take advantage of this informed new feature.
The remainder of this paper is organised as follows. Section 2 introduces the physical and mathematical background emphasising on signals projection modelling onto the coarray, as well as the stochastic behaviour of the envelope in tissue onto this same domain. Section 3 details the proposed beamformer: the Rician Beamformer (RiB). Section 4 describes the materials used and methods followed along the different tests described in Section 5. Finally, Section 6 presents the conclusions of this work.

Physical Model and SAFT-Based Coarray Projection
For the stochastic characterization of the underlying physical process, we start from Born model, notwithstanding its reported inaccuracies [15,19]. Different tissues are modelled as reflective point sources, grouped or not, in high density distributions, within an isotropic minimally dispersive medium. The structure is found in semi-space z > d min , whereas ROI is limited to plane y = 0-see Figure 1a. Regarding instrumentation, transducer has been modelled as a full multiple-input, multiple-output (MIMO) uniform linear array (ULA). It consists of N reciprocal sources, highly directive on φ and cuasi-omnidirectional on θ. Array elements are interleaved by λ/2. Electrical input to this array displays analytical form s(t) ∈ C, time separable and band-limited. For simplicity we consider a modulated Gaussian pulse with central frequency f 0 and bandwidth BW. Furthermore, the proposed methodology and conclusions are extensible whenever prior restrictions hold.
Hereafter, the bank of received signals is modelled as where F is the particle number, g i,j (t, p f ) is the result of the cascade filters modelling electroacoustic responses of the pairs, as well as the temporary space round trip of source f placed at p f . Considering Born model, (i, j)-th signal could be approximated as where c is the medium propagation speed and d x i ,p f ,x j is the round trip for transmitter i, reflective point source and receiver j-see Figure 1b. Reflectivity of a point source is represented by ρ f ∈ C. F ⊂ F is the subset of the sources close to y = 0 plane. On weakly dispersive media, it can be argued that fors(t) temporal separability is preserved, as g i,j (t, p f ) can be approximated by a delayed double derivative process. As SAFT is being used to compose the ultrasonic image, the pixel at p a underwent a focusing process at t = d x i ,p a ,x j /c on signals s i,j in S dataset.
To simplify Equation (2) UV coordinate system has been used. This introduces the concept of coarray [20]. After these transformations, the following approximation applies, spherical coordinates where the point source is placed and χ = x i + x j . The additive term ξ f accounts for errors due to far field approximation at point source f . As a result, the nearer the particle is with respect to the array, the larger the influence of ξ f .
In contrast to ULA, the χ codomain is unevenly sampled, whereas χ values are spatially sampled at λ/2. Both are related by a surjective function: where N is the number of elements of ULA, indexed by set N = {0, 1, . . . , N − 1}, considering a pulse-echo working mode, FMC signals are indexed by e = (i, j) tuples part of N 2 = N × N . The coarray domain is indexed by group L = N ⊕ N , where ⊕ denotes Kronecker sum. In turn, L is the result of concatenation: L := {L 0 ++ L 1 ++ . . . ++ L 2N−2 }, where L p denotes the group formed by the combination pairs of elements from N summing p. Its cardinality is a particularization of the dice rule. In this case, where Π[p] := 1 if 0 ≤ p < N; otherwise, 0. Considering the indexing imposed by the coarray, the signal corresponding to analysis point p a projected over the same domain is estimated as Based on (2), discrete estimates for the coarray in (4) and previous approximation d x i ,p a ,x j ≈ 2R a − u a χ + ξ a , the projection of the coarray can be modelled as: According to Equation (5), the projection onto the coarray domain for a focused point p a can be formulated as the summation of warped and shifted versions ofs(t). A little jitter effect for emitter-receiver pair positions must also be considered, represented by γ i,j , being the result of the biases introduced by ξ f and ξ a . Both α and β depend on the relative positions between p f and p a . These terms tend to 0 if p a ≈ p f .
Estimating statistics leads to a pernicious situation in a double way. On one hand, ULA-based systems do not contribute to diversity as they are supposed to do: when working on a pulse-echo set-up, all signals belonging to L p group are closely related, as they represent the same equivalent spatial frequency. On the other hand, acoustic noise due to low α f component biases the DAS estimator, s(χ) average.
From Equation (5), one may conclude the following.
• High correlation between different α f and β f suggests reflector clusters presence in tissue structures.
• Very low frequency components in s(χ) represent high reflectivity areas close to p a . Actually, only zero-frequency contributor is related to tissue response at that point.

SAFT Stochastic Behaviour
In recent years, statistical studies on radio frequency (RF) signals and ultrasonic image themselves have attracted increasing interest as they relate stochastic process models [21] with tissue physical parameters. One of the magnitudes attracting special interest is the envelope, which has been related to particle spatial organisation and reflectivity [22][23][24].
Initial approximations were carried directly on B-Scan signals envelope E [22], but nowadays, work is around specific image regions or global appearance [25,26].
Starting from (2) and its analytic form, knowing that s(t) is a modulated pulse: where F denotes all insonificated scatterers near y = 0. In the general case, E matches a homodyned K-distribution [27,28]. This is a variant of the Gamma distribution that yields a complex fitting process [29,30], but describes a large number of possible scenarios, i.e., flexible modelling of X and Y. However, these fitting processes lower parameter estimation accuracy as scatterers number increases [27]. Table 1 summaries distribution functions usually employed in ultrasound imaging to characterise envelop behaviour. Parameters could be interpreted as coherent signal power ε 2 , variance of scatterer strength σ 2 and α the scatterer clustering degree. Table 1. Distribution functions used to describe soft tissues behaviour. Overview of models in [29].

Name Probability Density Function Asymptotic Values
These models can be extrapolated for DAS conformed ultrasonic imaging, as for SAFT the image results from the absolute value of focused samples mean at p a . In analytic form The approximated expected value is estimated by the standard average out of N 2 samples. Once focused, the casuistic related to the X and Y variables modelling is dramatically reduced. This is because the same focusing process centres the zone under study, in spite of the integration inherent to the wave propagation process into the medium. Basically, X and Y reflect geometrical characteristics of the problem, displaying three different trends: • High particle count regions with low reflectivity discrepancies. In this case, both variables can be modelled as Gaussian noise with variance higher than the variables mean. The coarray projection in this zone shows random values of α f and β f in Equation (5). As a result, E can be modelled as a traditional random walk, following a Rayleigh distribution.
• High particle count regions with diverse reflectivity values, in geometrical groups that yield to energy concentration near the analysis point, i.e., concavities. Unfortunately, X and Y behaviour cannot be considered as monomodal: α f and β f get recurrent values due to the structure being analysed. In this case, the speckle noise is much higher in the image due to the rise in s(χ)| p a of low frequency waves. These aberrations cause artefacts due to the effects of secondary lobes and main lobe low directivity on the resulting image.
• High particle count regions with high reflectivity values and high spatial density. In this case α f and β f tend to 0; that is, there is high coherence between wavelets involved in summation (5). X and Y resemble non-zero mean Gaussian variables, and E follows a Rician distribution [31].

Rician Beamformer
The variability of the stochastic behaviour of the random variables involved in the focusing and projection processes of the coarray demands the integration of all three scenarios, prior to further analysis. For this, we exploit the space diversity principle.
The latter is commonly used in radar environments, not new to ultrasonic ULA systems, but up to now was mainly used for side lobes' suppression [32]. Sparse coarray sampling has proved valuable on modifying ULA directivity properties [33][34][35][36]. This is an usual strategy for resource minimisation [37,38]-i.e., coarray elements decimation-or to formulate effective acquisition strategies, at the expense of a minimum performance loss. All these approaches struggle between resolution improvement and lobes artefacts minimisation.
Space diversity is here exploited for sparse coarray random selection and applied on projected data, see Equation (5): First side lobe suppression criterion has been considered to determine required sparse level [39]. In this particular case, the projection s(χ)| p a is a 1.5 periods wavelet. A new random variable Z p is defined per each analysis point: By virtue of the central limit theorem, once a large enough number N z of virtual arrays groups A k have been formed, averages for the real and imaginary parts are equivalent to a pair of Gaussian random variables, identically distributed, with potentially non-zero mean.
This implies that in the three proposed scenarios in Section 2.2, random variable Z p can be modelled as a Rician, and fitting of Z p by means of the maximum likelihood method [40] provides an unbiased estimate of the Rician distribution. Figure 2 shows a naïve example of several image points for a one spot source case. For all evaluated cases measured contributors are due to the electrical and acoustic noises, except for the first case.

•
Case 1: Real target, hence signal to noise ratio (ε/σ) is high. In this scenario Z p acts as a Gaussian with √ ε 2 + σ 2 mean and σ 2 variance. In this specific case, µ| p a → ε.
• Cases 2 and 3: Most part of the noise is of acoustic nature and (ε/σ) ratio is low. Depending on the nature of α and β, see (5) (ε/σ) may be biased due to acoustic noise. Attending to (9), µ| p a → ε and σ 2 represents noise power and can be used as a penalty.

•
Case 4: Signal to noise ratio is low and most part of the noise contribution is due to thermal-electrical causes. These points can be used to determine the inherent noise power, using σ 2 estimator.
From these considerations, a new pixel-wise estimator is proposed for ultrasonic imaging, based on (10) fitting: where F Z p is the cumulative distribution function for Z p variable. The κ parameter controls penalty imposed by the proposed local expander-see Appendix A. Figure 2 depicts results from the proposed and DAS beamformers with κ = 0.01 for each of the last mentioned cases. For simplicity, Figure 2b depicts Z p fitting applied to the four considered cases, normalized by their maximum values. As the ultrasonic image is of logarithmic nature, to obtain a more convenient depiction, abscissa-axis is also logarithmic. Looking at the RiB(1%) results, the imposed penalty to low signal to noise ratio (SNR) points can be appreciated. One of the main advantages of the proposed method is local contrast ratio preservation in a pseudo-lineal way, while it also mitigates the speckle effect, especially in cavities and low reflection zones (Sections 5.3 and 5.4). The pseudo-lineal ratio between the analysis area and the resulting contrast depends on κ. By means of this ratio, beamformer expansion effect can be tuned, allowing its use on other applications such as NDT. As a counterpart, calibration process using standard phantoms is required, evaluating different κ values on testing cysts with fixed reflectivities, until converging to the expected phantom values, see Section 5.2.
Starting with the typical ultrasonic imaging model, where convolutive noise is present [11], adapted to TFM: where r DAS (x, z) is the resulting image analytic form, T is the Cartesian to prolate spheroidal coordinates transform, h(x, z) is the tissue response, n(ρ, ν) is the convolutive noise kernel and n a (x, y) is the set-up thermo-electric noise. The same aberration type holds in the resulting image using RiB, as long as κ 1: Nevertheless, one may argue that kernel memory for RiB,ñ(ρ, ν), is lower than n(ρ, ν), specially along ν dimension, alleviating speckle effect. On the other hand, additive noises, n a (x, z) andñ a (x, z), are seemingly the same.
This new ultrasonic imaging model is advantageous for most of despeckle techniques proposed on the last decade [12]. Despeckle tasks have to trade-off between analysis window size and image data loss [41]. A smaller acoustic noise kernel memory allows a smaller filtering window resulting in better image quality. This is specially noted on nonlinear and diffusion filtering.

Materials and Methods
We devised a number of tests to assess the proposed despeckle method. These we use to compare performance with state-of-art procedures, while covering a number of illustrative cases. In particular, the proposed experiences focus on three different questions: (i) axial and lateral resolution, (ii) correlation between data contrast obtained with the proposed method and the expected standard phantom values and (iii) reliability of soft tissue interfaces in complex structures.
Depending on the experiences and the specific setup, different arrays were chosen as listed in Table 2. Table 2. List of arrays employed in the following experiences. For real measures, acquisition equipment was a multichannel generation-acquisition system [42], developed by DASEL S.L. company (www.daselsistemas.es). A bandwidth of up to 68% has been achieved on all acquisitions. A sampling frequency of 40 MHz was set for all acquisitions using Array A, and 80 MHz otherwise-Array B.

Array Label Trademark Central Freq. (MHz) Number of Elements
For simulations (Array B) we used a 68% Gaussian pulsed signal and a white noise equivalent to 1 bit in a 12 bits sample depth. Array elements are 0.1386×1.386 mm in size and interleaved 0.1540 mm. For RiB current implementation, we chose N z = 10 4 virtual arrays.
Hereafter, we describe the devised tests, and compare both the proposed method (RiB) and DAS beamformer under different conditions.

High Reflectivity Threads Scenario
First, a simple test has been proposed in order to demonstrate the speckle effect reduction. This kind of test is related to the behaviour of point reflectors, but also can be used as a figure of merit to evaluate the procedure quality. To fulfil this task, a phantom has been made. It is composed of two nylon thread rows, normal to y = 0 plane, 0.35 mm width and interleaved 3 mm.
The scope of this experience is the evaluation of the RiB method on point sources (high reflectivity threads). RiB may yield to a better estimation of the interface impedance of the different tissues structures.

Standard Phantoms Contrast
The standard phantom model 040GSE (CIRS TM ) is a multipurpose multi-tissue phantom that includes several regions with distinctive reflectivity properties. This feature makes it ideal to test for any method to have a good contrast response. In this way, the proximal sets of grey cysts will be considered. This standard scale allows calibration, needed to compute κ actual value. Calibration procedure is common to all the samples analysed using the same set-up, as this threshold is mainly dependent on instrumentation features, due to system bandwidth, and is not so dependent on the sample itself.
As this experiment is not specifically designed to analyse the likelihood of the results compared to the expected values, a new study has been designed. This way several regions of the phantom have been analysed, following the manufacturer recommendations. Four different cysts have been considered with assigned reflectivities ranging from −6 dB to 6 dB with a diameter of 8 mm. Reflectivity values measured in the ultrasonic images are compared to background level on regions out of the cysts.

Tissue Interface Sharpness on Complex Phantoms
The proposed phantom in this section poses a challenge for ultrasonic imaging based on ULAs. The proposed simulation software for this task is Field II [43]. Array B has been used as transduction system, see Table 2. Two orthogonally interlocked hollow tori have been chosen as phantom geometry. Particle count across set-up being 10 6 , with uniform spatial density distribution.
Hollows of both tori are anechoic, whereas tori walls display reflecitivity 8 dB above background. This scenario is a particularly challenging one on the analysis plane (y = 0), as a rapidly changing, spiric section is found along this dimension. In Figure 3, a schematic representation of the proposed phantom is shown. Also, the theoretical image corresponding to y = 0 slice can be seen (projected over the y = −50 mm plane), as well as the slice plane itself, where the array will sense the phantom.  RiB and DAS results have been chained to non local means-based filtering (NLmeans) [44] as an example of RiB beamformer abilities and future potential, although other more complex non-linear filters can be tested [45]. Filter parameters tuning are 5 for search window, 3 for similarity window and a filtering order of 10.

Stationary Heart Standard Phantom
Finally, one of the most popular applications of ultrasonic imaging segmentation is heart cavities study [46]. This problem is specially sensitive to speckle noise, due to the large impedance differences between walls and hollows, causing considerable segmentation errors.
In the aforementioned tests, RiB effectiveness has been demonstrated, via simulation, on hollowed complex form impedance analysis. In this section, the study will be extended to real cases. To do that, Ultrasound Heart Phantom model 067, CIRS TM , will be employed, jointly to the same experiment set-up using Array A (see Table 2). This phantom has been designed for cardiac imaging training.
As in the last test, the proposed ultrasonic images have been calculated-DAS y RiB(10%)-as well as NLmeans filtering results, using the same parameters as in Section 5.3.

Results
Hereafter, we present results obtained both from real measurements and simulations to prove the robustness of RiB method.

High Reflectivity Threads Scenario
For this test, Array A has been used (see Table 2), tuned as indicated in the table. In Figure 4   Using the proposed method (RiB 10% and RiB 1%, Figure 4b,c), a more reliable image with respect to the analysed structure is obtained, with a considerable increase on image contrast. A point to point comparison of the values obtained using DAS vs RiB methods is depicted in Figure 4f. Theoretical 1:1 ratio is included in red. As illustrated in Figure 4a-c the more reflective the zone, e.g., upper part of graph in Figure 4f, all beamformers behave similarly. However, on regions corresponding to the transmission medium, these suffer a considerable attenuation, inversely proportional to κ.

Standard Phantoms Contrast
The ultrasound image of the phantom is shown in Figure 5. Cysts presence as long as the expected reflectivity can be observed in both results-DAS: Figure 5a and RiB: Figure 5b. However, the contrast of the proposed beamformer is far better.
The primary objective of the experiments is to determine the most accurate value for κ to achieve the best ultrasonic image fitting to the expected reflectivity values provided by the phantom manufacturer. As the image resulting from RiB displays linear behaviour with respect to log κ, see Appendix A, the tested values follow the same progression. The selected value for RiB image processing is κ =10%. Based on this classification and the parameter selected value, the result of the analysis for both beamformers is shown in Figure 5c in form of a boxplot. This way outliers for each region are also displayed. Boxplots are paired, the leftmost is the result of DAS, whereas RiB output are the rightmost. Focusing on the difference of the median values for the regions, compared to the median value of the background, contrast ratio gain is evidenced. Several conclusions can be obtained: both methods display progressive scale of grey values corresponding to each reflectivity expected value for the cysts. However, DAS beamformer displays no sensible difference while comparing 3 dB and 6 dB regions. Difference between regions measured using RiB beamformer displays a 3 dB variation, in coherence with the manufacturer specifications. The interquartile range (IQR) of the RiB regions is wider because this method enhances texture displayed and contrast.

Tissue Interface Sharpness on Complex Phantoms
Images in Figure 6a,b correspond to the results of DAS and RiB on the simulated phantom, respectively. The most significant difference between both is the sharpness of the spiric section. This is mainly due to the lowering influence of the tissue interface close to y = 0 plane. The speckle introduced by that interface is considerably reduced after RiB processing. Also, due to the overall image contrast improvement, anechoic hollows are clearly distinguished. As in the last results, RiB shows better performance regarding tissue interface sharpness and anechoic regions characterisation, see Figure 6d. In fact, after NLmeans filtering, RiB result poses a successful approximation to the expected spiric section. Contrarily, filtering result on DAS, see Figure 6c, flattens the spiric section and the deeper parts of the torus in the centre of the image. This is clearly rejected in the case of RiB. This brings to light differences on kernel memories for convolutive noise-see (12) and (13).

Stationary Heart Standard Phantom
Finally, the RiB method was tested on an actual measurement with notable impedance differences. One of the most obvious scenarios for this kind of test is cardiac imaging. In this case, the same before mentioned setup has been used, resulting on images in Figure 7. Comparing DAS, see Figure 7a, to RiB, see Figure 7b, one may observe how speckle effect has been considerably mitigated. RiB speckle reduction is especially evident on cavities where DAS reflects reverberations-convex interfaces with respect to the array.
Conversely, in Figure 7c,d NLmeans filtered image versions are depicted, using the same parameters as in Section 5.3. For these parameters setup and the resolution used to generate the images, search window is approximately 1.6 mm. This, and filter memory used to model convolutive noise, see (12) and (13), explain aberrations on DAS filtered image, Figure 7c, hindering heart cavity segmentation. RiB definition facilitates correct segmentation, delivering good tissue interface discrimination capabilities, a desirable feature for this kind of medical imaging applications.

Discussion
One of the most disturbing phenomena affecting ultrasonic image quality is nonuniformity of noise variance along the inspection area. This compromises CAD tasks and image readability. The proposed method mitigates this phenomenon from two different perspectives:

•
Contrast improvement. RiB improves image contrast without affecting texture. Compared with standard phantom data, contrast values obtained from RiB are consistent.
• Tissue interface sharpness. RiB notably improves tissue interface contrast, especially in the case of tissue-cavity interfaces, in particle-free cavities.
Results obtained by means of RiB beamforming can be fed to further denoising systems, highly improving results. In this paper, pdf usage has been restricted to a single point evaluation related to DAS value. Point-wise, RiB estimations could be combined with alternative sources of information to produce enhanced images (fusion), easier to interpret.
Robust mean and variance estimations are obtained with this method, thanks to diversity principle, applied to coarray projection. Estimations can be used to implement new image domain despeckle methods that exploit noise variance variability. Such information is instrumental to more complex denoise filtering methods, and also could be used in other medical imaging areas, such as Magnetic Resonance Imaging (MRI). Though it is out of the scope of this contribution to formulate this kind of filtering processes, there is increasing interest among the scientific community in applying these methods to tissue medical imaging.

Conflicts of Interest:
The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. RiB First-Order Cumulative Distribution Functions
Let X be a random process for which first-order probability distribution matches Rician distribution with parameters ε and σ, We define ρ = ε/σ, essentially being a signal-to-noise ratio, and introduce the corresponding cumulative distribution function, defined as On this general form, one may look into particular cases described in terms of distribution parameters, and interpret the RiB definition in (4). Penalties corresponding to the Rician distribution are depicted in Figure A1. This is, the difference between DAS and RiB outputs. Parameter κ significantly affects low SNR values, while it preserves high SNR inputs. Whenever ε → 0, this implies that ρ → 0 and I 0 xε σ 2 → 1 for any possible value of x. In this case, no coherent contributors arise in (5).