3D Synthetic Aperture Imaging Method in Spectrum Domain for Low-Cost Portable Ultrasound Systems

Portable, hand-held ultrasound devices capable of 3D imaging in real time are the next generation of the medical imaging apparatus adapted not only for professional medical stuﬀ but for a wide group of less advanced users. Limited power supply capacity and the relatively small number of channels used for the ultrasound data acquisition are the most important limitations that should be taken into account when designing such devices and when developing the corresponding image reconstruction algorithms. The aim of this study was to develop a new 3D ultrasound imaging method which would take into account the above-mentioned features of the new generation of ultrasonic devices – low-cost portable general access scanners. Itwas based on the synthetic transmit aperture (STA) method combined with the Fourier spectrum domain (SD) acoustic data processing. The STA using a limited number of elements in transmit and receive modes for ultrasound data acquisition allowed both aforementioned constraints to be obeyed simultaneously. Moreover, the computational speed of the fast Fourier transform (FFT) algorithm utilized for the ultrasound image synthesis in the spectrum domain makes the proposed method to be more competitive compared to the conventional time domain (TD) STA method based on the delay-and-sum (DAS) technique, especially in the case of 3D imaging in real time mode. Performance of the proposed method was veriﬁed using numerical 3D acoustic data simulated in the Field II program for MATLAB and using experimental data from the custom design 3D scattering phantom collected by means of the Verasonics Vantage 256™ research ultrasound system equipped with the dedicated 1024-element 2D matrix transducer. The


Introduction
Medical ultrasonography is the most used imaging modality for its speed, flexibility, cost-effectiveness and non-invasive nature.These features make the ultrasound imaging more competitive than the other imaging techniques, especially for preliminary diagnosis.Recently, the 3D or volumetric ultrasound imaging technology has been rapidly developed as the researchers started exploring innovative and new applications (Fenster et al., 2001;Ji et al., 2011).The volumetric imaging provides the physicians with much deeper and realistic insight into the examined part of the human body (Campbell et al., 2005;Landry et al., 2005).For instance, 3D ultrasound provides the capabilities of surface characterization of the fetus offering a better view of fetal defects like a cleft lip or clubfoot.Also, reduction in breast examination time with 3D ultrasound is currently under investigation (Kotsianos-Hermle et al., 2009; Padilla et al., 2013).
The most recent advances in ultrasound technology are enabling the rapid development of 3D imaging in real-time allowing for volumetric visualization and tracking the movements of imaged organs, like the motion of the human heart wall or valves, blood flows in various vessels and so on.Moreover, the real-time 3D ultrasound can be potentially integrated into different portable small-size devices suitable for diagnostic and therapeutic procedures (due to small probe size and safety) in the near future.
There are a number of constraints that should be considered when designing such devises, like the limited power consumption due to the limited capacity of the power supply battery or a huge amount of acoustic data that have to be acquired using a limited number of signal channels and processed on-the-fly during volumetric data acquisition with 2D matrix transducers consisting of thousands of elements and so on.
To mitigate this problem the sub-array architectures were considered in the literature based on the sub-aperture design (Savord, Solomon, 2003) or a separable beamforming design wherein a conventional 2D array beamforming is decomposed into a series of 1D beamforming problems at the cost of some reduction in image quality (Yang et al., 2013).Another strategy to reduce the number of active elements is by using the sparse 2D arrays.For example, Austeng and Holm (2002) proposed a method based on the principle of grating lobes suppression to form the optimal symmetric and non-symmetric regular sparse periodic and radially periodic designs.Sparsely sampled aperture patterns like the minimalredundant arrays were proposed in (Karaman et al., 2009).To reduce redundancy the authors developed a sparse array design procedure based on the spatial convolution of the transmit and receive patterns by using different combinations of linear sub-arrays of the 2D matrix transducer.Several non-rectangular 2D sparse array layouts have also been proposed, such as the concentric circular arrays (Wang et al., 2002;Ullate et al., 2006) and spiral arrays (Martínez-Graullera et al., 2010; Yoon, Song, 2019).A more detailed literature review on the 2D sparse arrays design can be found for example in (Ramalli et al., 2022).
Recently, the spectrum domain (SD) ultrasound imaging methods have been attracting growing attention due to their high computational speed.This is a significant advantage over the conventional time domain (TD) methods based on the classical delay-andsum (DAS) technique (Thomenius, 1996).This makes them promising for the 3D ultrasound imaging in the real-time mode.Several SD methods have been reported in the literature (Busse, 1992;Skjelvareid, 2012;Cheng, Lu, 2006).In (Cheng, Lu, 2006), theoretical development of the so-called phase migration method was presented.The data collected using 2D aperture were first transformed for each axial depth to the spectrum domain.Then the frequency-dependent propagator function was applied to each cross-sectional spectrum followed by the frequency averaging within the considered frequency band.The magnitude of the inverse Fourier transform of the averaged spectra obtained for each axial depth independently, yielded the 3D ultrasound image.The numerical complexity of the phase migration method was comparable to the conventional 3D TD imaging methods at the same time this yielded the synthesized images of poor quality, therefore remaining a purely theoretical interest.In (Skjelvareid, 2012), the well-known ω − k migration algorithm previously adopted for the 2D B-mode ultrasound (Skjelvareid et al., 2011) was generalized for the 3D imaging.The ultrasound data were first transformed to the spectrum domain.Next, the ω − k migration was applied to transform the temporal spectral variable ω to the spatial spectrum variable related to the axial spatial coordinate.Then the inverse 3D Fourier transform was computed to obtain the final 3D ultrasound image.In (Cheng, Lu, 2006) yet another approach for 3D SD ultrasound image synthesis was proposed.It was based on the assumption of the plane acoustic wave insonification (Montaldo et al., 2009) by an infinite 2D aperture array.The acoustic signals detected by the same 2D array were represented as the superposition of harmonic waves reflected from the scatterers distribution in the examined volume.Each harmonic component within the limited frequency band was first transformed to the spectrum domain with respect to the lateral coordinates.Next, the mapping between the wave-vector components of the detected echoes and the back-scattered field related to the reflectivity function describing the scattering properties of the examined medium was applied in the spectrum domain.Finally, the inverse 3D Fourier transform yielded the final 3D ultrasound image.
The main objective of this research was to propose a new spectrum domain synthetic transmit aperture (SD-STA) 3D image reconstruction method.It offers a promising solution to the aforementioned constraints, suitable for low-cost portable devices which should operate with a limited number of signal channels and limited power consumption.Specifically, the SD-STA method is based on a novel approach which combines the non-overlapping sub-aperture data acquisition of the back-scattered ultrasound echoes by analogy with the TD-STA methods (Gammelmark et al., 2003;Jensen et al., 2006;Nikolov et al., 2008;Trots et al., 2009;Tasinkevych et al., 2012) and the ultrasound data processing in the Fourier spectrum domain.Using a limited number of elements in the transmit (TX) and receive (RX) modes allowed the amount of acoustic data collected and transferred during a single TX/RX event to be reduced by utilizing a limited number of signal channels (64 and 256 considered in this study).Moreover, the computational speed of the fast Fourier transform (FFT) algorithm used for the ultrasound image synthesis in the spectrum domain makes the SD-STA method promising for the real-time ultrasound imaging.
The validation of the proposed method was carried out in two stages.First, the method was tested using numerical acoustic data simulated in MATLAB using Field II software (Jensen, 1996;Jensen, Svendsen, 1992).Next, the method was validated using experimental acoustic data acquired in a custom made 3D scattering phantom by means of the Verasonics Vantage 256™ research ultrasound system equipped with the dedicated 1024-element 2D array transducer operating at 3.47 MHz center frequency.The comparative analysis of the proposed SD-STA method with the conventional TD-STA method (Tasinkevych et al., 2013) generalized for the case of 3D imaging (Tasinkevych, 2017), was conducted.It was evidenced both in simulations and measurements that the proposed SD-STA method provides high computational speed and imaging quality comparable to stateof-the-art TD-STA methods.
The rest of the paper is organized as follows.In the next section theoretical background of the proposed SD-STA method is given.Then, the method is described in Sec. 3. The results of 3D imaging using numerically simulated and experimentally obtained data are presented in Sec. 4. The discussion of the results obtained is given in Sec. 5. Finally, in Sec.6 the summary of the work is presented.

Theory
Theoretical backgrounds of the SD-STA methods are briefly presented in what follows.Let one consider an N × N -element flat matrix transducer placed in the plane z = 0.The origin of the Cartesian coordinate system is in the center of the transducer aperture and its sides are parallel to the x-and y-axes.During a single TX/RX event the i-th M × M -element sub-aperture is excited with a short pulse (one sine cycle of the nominal frequency) and transmits an unfocused wave in the z-axis direction, where i = 1, ..., I, I = (N M ) 2 .Without loss of generality the back-scattered echoes can be assumed to be detected by the same M × M -element sub-aperture and the i-th LRI is then synthesized using the RF echoes detected during this single TX/RX event only (this case is referred to as the reduced matrix of M × M RF echoes in the Subsec.3.1).A more general case of using several TX/RX events for the synthesis of a single LRI is discussed in the next section (this case is referred to as the full-size matrix of N × N RF echoes in Subsec.3.1).The received echoes represent the spatially sampled back-scattered acoustic field s (r , t) in the plane z = 0 within the 2D sub-aperture; r ≡ (e x x + e y y + e z 0) is a position vector defining a location of individual element of the transducer in the plane z = 0; s (r, t) is the back-scattered acoustic field in the volume of interest V ; r ≡ (e x x + e y y + e z z) is a position vector in V .The spatial spectrum of the back-scattered filed detected by the i-th RX subaperture can be represented as a linear combination of harmonic components propagating form randomly distributed point scatterers in the volume V (Cheng, Lu, 2006): where f (r) denotes an object function describing the scattering intensity of the random distribution of scatterers in the volume V .In Eq. ( 1) x + e y k I y + e z k I z is the wave-vector of the transmitted acoustic wave.In the case of SD-STA considered in this paper k I = e z k I z ; also, in Eq. (1) k = e x k x + e y k y + e z k z is the wave-vector of the backscattered wave; k = e x k I x +e y k I y and K(ω) is a transfer function that accounts for the influence of the transmit and receive circuitry.From Eq. (1) the relationship between the temporal-spatial spectrum of the backscattered acoustic field S (k , ω) and the spatial spectrum of the object function F (k ′ ) in V can be obtained: where F ′ (k ′ ) i is the so-called band-limited spatial spectrum of f (r) obtained as a result of transmission and detection of the acoustic wave by i-th M × M -element sub-aperture.It defines the 3D spatial spectrum of the partial LRI reconstructed using the RF echoes detected by the i-th M × M -element sub-aperture: The spatial spectrum defined in Eq. ( 3) can be obtained from the temporal-spatial spectrum of the detected back-scattered echoes as: where the spatial spectrum S (k ′ (ω)) i of the backscattered signals is obtained from S (k , ω) i after the variable transform: To obtain the final 3D HRI of the examined volume of interest its 3D spatial spectrum has to be synthesized first as a sum of the spatial spectra of all LRIs: followed by the inverse 3D Fourier transform applied to I (k ′ (ω)): In Eq. ( 7) the term B (k ′ (ω)) was introduced which is the Jacobian of the variable transformation in the spectrum domain defined in Eq. ( 5):

Methods
The SD-STA method proposed in this paper was tested using numerical data simulated in MATLAB using Field II software.Moreover, experimental measurements were conducted and acoustic data from the custom design 3D scattering phantom were collected using Verasonics Vantage™ research ultrasound system equipped with a 2D matrix transducer utilizing 1024 elements in a 32 × 32 grid and operating at 3.47 MHz center frequency.

Numerical simulation in Field II
First, the method was tested using numerical data simulated in Field II.For this purpose the 32 × 32element matrix array transducer operating at f 0 = 3.47 MHz frequency was modeled.The transducer element pitch was 0.3 mm and the element width was 0.275 mm.The sampling frequency was 13.87 MHz which corresponded to 4 time samples per acoustic wavelength of the recorded signal.There were N t = 1122 time samples simulated in each RF echo (for the 60 mm depth assumed, see further discussion).The speed of sound c = 1540 m/s and the frequency dependent attenuation α = 0.5 dB/[MHz ⋅ cm] around the center frequency f 0 (Jensen, 2021) were applied which is typical for the soft tissue (Hill et al., 2004).The parameters applied in the numerical simulations were chosen to mimic the 2D matrix transducer used in the experimental measurements (see discussion in the next section).The transducer was excited with a short pulse burst (one cycle of the nominal frequency).The sets of 5 anechoic spheres (cysts) and 5 scattering spheres located in the volume z × x × y = 10 × 10 × 60 mm filled with randomly distributed point-like scatters were simulated separately.The diameter of spheres was 5 mm.The spheres were spaced 10 mm axially and located on the z-axis normal to the 2D aperture.The number of point scatterers was 30 per mm 3 (in accordance with the Rayleigh scattering conditions), yielding the total number of 1.8e5 point-like scatterers in the volume under consideration.In the case of cysts the uniform scattering amplitude was assumed for the scatterers distributed outside the cysts and the scattering amplitude for the scatterers inside the cysts was set to zero (anechoic spheres).In the case of scattering spheres the ratio of the scattering amplitudes for the inner and outer scatterers was assumed 10:1.Numerical simulations were conducted for nonoverlapping TX/RX sub-apertures arranged in 16 × 16 and 8 × 8-element grids.Two different data acquisition schemes were examined.They differed in the number of RF echoes collected for a given TX sub-aperture which were then used for a single LRI synthesis.
First, to synthesize a single LRI corresponding to a given TX sub-aperture the RF echoes were simulated for all elements of the 2D matrix transducer, as illustrated in Fig. 1 for the particular case of 16 × 16element TX/RX sub-aperture.
Specifically, to synthesize the LRI #1 the interrogating pulse was transmitted four times by the TX sub-aperture #1 for each of the four RX sub-apertures used successively in the RX mode which yielded the full-size matrix of 32 × 32 RF signals denoted by [RX1] in Fig. 1.This procedure was then repeated for each of the four TX sub-apertures.All the LRIs synthesized using the full-size matrices [RX1] through [RX4] of the RF signals collected this way were then summed in the spectrum domain yielding the final HRI.Therefore, to obtain the final image there were 16 TX/RX events required for the case of the 16 × 16-element sub-aperture data acquisition and 64 TX/RX events for the case  of the 8 × 8-element one, respectively.These cases were denoted as TX16/RX32 and TX8/RX32 in Sec. 4, correspondingly.
Next, to obtain a single LRI corresponding to a given TX sub-aperture the so-called reduced matrices of the RF signals were simulated.In particular, for the case of 16 × 16-element TX/RX sub-apertures to synthesize the LRI #1 corresponding to the TX subaperture #1 the interrogating pulse was transmitted only once and the RF echoes were then collected using the same grid of elements switched to the RX mode as depicted in Fig. 2.This yielded a reduced matrix [RX1] of 16 × 16 RF signals (see Fig. 2).
The rest of the full-size 32 × 32 matrix corresponding to inactive elements of the 2D array transducer was filled with zeros (as discussed later in this section) and the LRI #1 was synthesized then in the spectrum domain.This procedure was repeated for each of the For j-th M × M -element TX sub-aperture, j = 1, ..., J, J = (N M ) 2 : Reduced RF matrix Full-size RF matrix 1) Transmit interrogating pulse with j-th M × Melement TX sub-aperture and acquire ultrasound echoes using the j-th M × M -element RX sub-aperture.
2) Fill the rest of N × N RF matrix with zeros (empty RFs).
1) Transmit interrogating pulse with j-th M × Melement TX sub-aperture and acquire ultrasound echoes using the 2) Fill full N × N RF matrix with acquired echoes.
3) Compute the 3D temporal-spatial spectrum S (k , ω) j of the echoes with respect to the spatial variables x, y and the time t using the FFT algorithm.4) and ( 5)) to obtain the spatial spectrum I (k ′ ) j of the j-th LRI. 5) Accumulate the spatial spectrum of the HRI: I (k ′ ) = I (k ′ ) + I (k ′ ) j (see Eq. ( 6)).

4) Transform the variables
four TX sub-apertures and the final HRI was obtained in the same way as in the case of the full-size RF matrix data acquisition discussed earlier.Hence, to obtain the HRI in the case of the reduced RF matrix approach shown in Fig. 2 there were only 4 TX/RX events required for the 16 × 16-element TX/RX subaperture data acquisition and 16 TX/RX events for the 8 × 8-element one, respectively.These cases were denoted as TX16/RX16 and TX8/RX8 in Sec. 4, correspondingly.
Finally, for comparison the simulations were also conducted for the full-size 32 × 32-element TX/RX aperture.In this case a single interrogating pulse was transmitted and the full-size matrix of 32 × 32 RF echoes was collected at once.These echoes were then used to synthesize the final HRI in the spectrum domain.This case was denoted as TX32/RX32 in Sec. 4.
The data-flow of a single LRI synthesis is sketched: The ultrasound data needed to synthesize a single LRI were acquired first yielding the full-size or reduced matrix of RF echoes (step 1 and 2).Specifically, in the former case the data acquisition for the j-th M × M -element TX sub-aperture yielded a 3D matrix of N × N × N t time samples, where N t is the number of time samples in a single RF echo (N t = 1122), and N × N is the 2D matrix transducer size (N = 32).In the case of the reduced RF matrix approach only an M × M × N t sub-matrix of time samples was acquired during the j-th TX/RX event.The rest of the transducer elements remained inactive.To synthesize the partial LRI in this case the rest of the 3D matrix of N × N × N t time samples was filled with zeros which corresponded to empty RF echoes (N t zero samples in each of them).This allowed the spectrum domain data processing to be more versatile and suitable for different TX/RX sub-aperture size and position within the 2D transducer array which greatly simplified the HRI spatial spectrum synthesis (step 5).The 3D matrix of time samples was then used to compute the spatialtemporal spectrum S (k bot , ω) with respect to the lateral spatial variables x and y and the time variable t (step 3).For this purpose the fft.m routine from the MATLAB Signal Processing Toolbox was used.The variable transformation in the spectrum domain in the step 4 was implemented as a 1D interpolation procedure.Concretely, the temporal spectrum S (∶, ω) j was interpolated to the spatial spectrum I (∶, k z (ω)) j using Eq.(5).For this purpose the routine interp1.m from MATLAB was used.In the final step 5 the spatial spec- trum I (k ′ ) j of the j-th LRI was added (accumulated) to the spatial spectrum I (k ′ ) of the final HRI.Once all I (k ′ ) j were synthesized the inverse 3D Fourier transform of I (k ′ ) yielded the final 3D ultrasound image I (r) (see Eq. ( 7)). Figure 3 depicts the comparison of ultrasound image formation for the SD and TD methods considered in this paper.
It should be noted that in the case of the TD-STA method the synthesized HRI is represented by the properly delayed and summed RF echoes collected during consecutive TX/RX events (Tasinkevych et al., 2013).Therefore, the envelopes of the synthesized RFs must be computed first as shown in the diagram in Fig. 3.For this purpose the routine envelope.mfrom the MATLAB Signal Processing Toolbox was used.Furthermore, for consistency the images obtained using the TD-STA method were displayed as a function of spatial variables as in the case of the SD-STA one.For this purpose the relation between the axial spatial variable z and the time t was used: z = ct 2, c being the acoustic wave speed.
In the numerical simulations the 3D ultrasound images comprised of N x × N y × N z = 64 ×64 ×512 spatial samples were synthesized using the SD-STA method proposed.Also, the same acoustic data were processed using the TD-STA method (Tasinkevych et al., 2013; 2017) and the corresponding results were compared.
The cysts arrangement considered in the numerical simulations enabled the comparison of the image contrast and contrast-to-noise ratio (CNR) obtained for different data acquisition schemes considered.The contrast and CNR were defined as: where µ c , µ r are the mean intensities of cyst and reference regions; σ c , σ r are the standard deviations for the cyst and reference regions, respectively.The corresponding values of the contrast and CNR were assessed from the 2D B-mode images of the cross-sectional views of the cysts arrangement (the plane Y = 0).The cyst region in Eq. ( 9) was defined as a circle area of the diameter D c = 0.75D centered within the corresponding cross-section of the cyst, whereas the reference region was defined as an area between the circles of diameters 1.5D and 1.25D surrounding the cyst crosssection.Finally, the 3D imaging capability of the proposed method was demonstrated using acoustic data simulated for the set of scattering spheres arrangement discussed above.The isosurfaces at the fixed intensity level of the corresponding 3D reconstructed spheres were obtained using the sliceomatic visualization package for MATLAB.

Experimental measurements
Next, the SD-STA method was verified using measurement data.The experimental setup is schematically in Fig. 4. The ultrasound data for the scattering sphere with a diameter of 5 mm (shown on the right in Fig. 4) were acquired using the Verasonics Vantage 256™ research system (Kirkland, WA) equipped with the dedicated matrix transducer utilizing 1024 elements arranged in a 32 × 32 grid and operating at 3.47 MHz center frequency.The examined sphere had a through hole with a diameter of about 1 mm.It was suspended in the water tank using a thin nylon wire with a diameter of 0.2 mm at a distance of 25 mm from the 2D array transducer face.The speed of sound c = 1490 m/s in water at the room temperature was applied in the image reconstruction (Hill et al., 2004).The excitation voltage was 10 V with 12 dB amplification of the received signals.The scanner was preprogramed to scan the phantom according to the 5 data acquisition schemes described in Subsec.3.1.It should be noted that in the case of the full-size 32 × 32-element TX/RX aperture the acoustic data were acquired using the dedicated software provided by the Verasonics Vantage 256™ research system manufacturer.It allowed a single flash transmit over all 1024 elements to be obtained, with each of the four 256-element RX sub-apertures used in the RX mode.Therefore, to collect the fullsize matrix of 32 × 32 back-scattered RF echoes in this case four consecutive TX/RX events were required.

Numerical simulations
The B-mode images of the numerical acoustics data simulated in Field II for the set of cysts (anechoic spheres) of 5 mm in diameter arranged axially in front of the 2D matrix transducer are shown in Figs. 5 and 6.The images in Figs. 5 and 6 present the crosssectional views of the 3D reconstructed data obtained using the SD-STA method proposed in this work (Fig. 5) and the conventional TD-STA method generalized for 3D imaging (Fig. 6).All B-mode images were In Figs. 7 and 8 the plots of the CNR and contrast assessed from the B-mode images depicted in Figs. 5  and 6, are shown, respectively.In Fig. 9 the 3D images of the numerical acoustics data simulated in Field II for the set of scattering spheres arranged axially in front of 2D matrix transducer are shown (see Subsec.3.1).The results were obtained using the SD-STA method for different data acquisition schemes.The images in Fig. 9 present the isosurfaces of the normalized 3D reconstructed data in logarithmic scale generated using sliceomatic visualization package for MATLAB and their arrangement is similar to that in Fig. 5.

Measurements
In Fig. 10 the 3D images of the experimental acoustics data acquired from the custom design scattering phantom (see Subsec.3.2) are shown.The results were obtained using the SD-STA method for different data acquisition schemes.The images present the isosurfaces of the normalized 3D reconstructed data in logarithmic scale generated using sliceomatic visualization package for Matlab.They are arranged in the same order as in Fig. 5.The B-mode images of the cross-sectional views corresponding to the X = 0 and Y = 0 planes of the experimental acoustics data acquired from the custom design scattering phantom are presented in Figs.11  and 12, respectively.All B-mode images were normalized with respect to their maximum value in the corresponding cross-sections and visualized in logarithmic scale over 30 dB dynamic range.
Fig. 11.B-mode images of the experimental acoustic data corresponding to the crosssectional views (plane X = 0) of the 3D reconstruction obtained using the FD-STA method for different data acquisition schemes (see Fig. 10).The images are plotted in logarithmic scale over 30 dB dynamic range.
Fig. 12. B-mode images of the experimental acoustic data corresponding to the crosssectional view (plane Y = 0) of the 3D reconstruction obtained using the FD-STA method for different data acquisition schemes (see Fig. 10).The images are plotted in logarithmic scale over 30 dB dynamic range.

Discussion
The results presented in this work confirmed effectiveness of the SD-STA 3D imaging method.Different data acquisition schemes based on the nonoverlapping sub-apertures used in the TX/RX modes were studied and compared.For this purpose the numerical acoustic data simulated in Field II program and experimental data collected with the Verasonics Vantage 256™ research ultrasound system from the custom-design 3D scattering phantom were used.It was evidenced in this research both by numerical simulations and experimental data imaging that the tradeoff exists between the frame-rate defined primarily by the number of acoustic data frames (volumes) acquired, and the imaging quality, assessed by the contrast and the CNR.Specifically, the B-mode crosssectional views (see Fig. 5 in Subsec.4.1) of the 3D numerical acoustic data simulated in Field II program for different data acquisition schemes and reconstructed using the SD-STA method were analyzed and compared.Only an insignificant change in the contrast and CNR for the depth below 40 mm was observed, provided the full-size matrices of 32 × 32 RF echoes was used for the synthesis of partial LRIs as evidenced further in this section.These cases were de-noted as TX16/RX32, TX8/RX32, and TX32/RX32 in Figs.7a and 8a, respectively.At the larger depths a slight increase in the contrast and CNR was observed for the TX16/RX32 and TX8/RX32 data acquisition schemes in comparison to the TX32/RX32 one.Specifically, about 9 dB and 5 dB increase in the contrast and 10 dB and 6 dB increase in the CNR was obtained for the TX16/RX32 and TX8/RX32 cases over the TX32/RX32 case, respectively.In the case of the TX16/RX32 data acquisition scheme, the final HRI was obtained by summing 4 LRIs which required 16 TX/RX events, while in the case of the TX8/RX32 one 16 LRIs and 64 TX/RX events were needed.Obviously, the time needed to acquire the acoustic data and synthesize the HRI determines the frame rate of the imaging method which is one of the most important parameters especially in the real-time 3D ultrasonography.Therefore, to reduce the number of TX/RX events required for the synthesis of the HRI the data acquisition schemes TX16/RX16 and TX8/RX8 were implemented and tested (see Subsec.3.1).Specifically, a single TX/RX event was needed to acquire the reduced matrices of 16 × 16 and 8 × 8 RF signals with non-overlapping TX/RX sub-apertures comprised of 16 × 16 and 8 × 8 elements, respectively.These reduced matrices were then used for the LRIs synthesis.The number of the LRIs needed to obtain the HRI did not change, but the data acquisition time was shortened 4 times for the TX16/RX16 and 16 times for the TX8/RX8 data acquisition schemes in comparison to the TX16/RX32 and TX8/RX32 ones, respectively.However, this increase in the data acquisition rate was achieved by the cost of the synthesized image quality deterioration, as evidenced by examples shown in Fig. 5, panels 4 and 5 denoted as the TX16/RX16 and TX8/RX8, respectively.This degradation of the image quality was quantified by the contrast and CNR, shown in Figs.7a and 8a.Specifically, the decrease in the contrast was 3-5 dB at the depths below 40 mm up to 14 dB at the maximum simulation depth of 50 mm and the decrease in the CNR was 4-7 dB at the depths below 40 mm up to 18 dB at the 50 mm depth for the TX16/RX16 data acquisition scheme in comparison to the TX16/RX32 one, respectively.The corresponding values of the contrast and CNR decrease for the case of the TX8/RX8 data acquisition scheme in comparison to the TX8/RX32 one were 3-7 dB (contrast, Fig. 8a) and 6-8 dB (CNR, Fig. 7a) at the depths below 40 mm and up to 13 dB (contrast, Fig. 8a) and 14 dB (CNR, Fig. 7a) at the depth of 50 mm, respectively.
Similar dependence of the image quality on the number of TX/RX events based on the visual assessment was also observed in the case of the experimental data reconstruction.The corresponding examples are shown in Figs.11 and 12, where the B-mode crosssectional views in the planes X = 0 and Y = 0 of the 3D reconstructed data of the custom design scattering phantom (see Fig. 10) are depicted.
Furthermore, the performance of the SD-STA method was compared with the conventional TD-STA ultrasound imaging method.The same TX/RX data acquisition schemes were implemented and the image quality was assessed for both methods using the contrast and CNR parameters.In the case of the TD-STA method similar dependencies of the contrast and CNR parameters versus depth for different TX/RX data acquisition schemes were obtained as in the case of the SD-STA method.For instance, for the TX16/RX16 data acquisition scheme the decrease in the contrast was 2-8 dB at the depths below 40 mm up to 20 dB at the maximum simulation depth of 50 mm in comparison to the TX16/RX32 one, respectively (see Fig. 8b).The corresponding decrease in the CNR was 2-11 dB at the depths below 40 mm up to 22 dB at the 50 mm depth (see Fig. 7b).In the case of the TX8/RX8 data acquisition scheme the observed decrease in the contrast and CNR in comparison to the TX8/RX32 one were 1-5 dB (contrast, Fig. 8b) and 1-7 dB (CNR, Fig. 7b) at the depths below 40 mm and up to 11 dB (contrast, Fig. 9b) and 9 dB (CNR, Fig. 8b) at the depth of 50 mm, respectively.
Although the qualitative dependencies of the contrast and CNR vs. depth for the SD-STA and TD-STA methods were similar, the quantitative changes in the values of corresponding parameters differed slightly for the data acquisition schemes considered.For example, comparison of the images reconstructed using both methods from the full-size matrices of 32×32 RF echoes confirmed that for the TX16/RX32 data acquisition scheme the contrast decreased approximately by 7 dB (from −2 to −9 dB, Fig. 8a) in the case of the SD-STA method and by approximately 4 dB (from −3 to −7 dB, Fig. 8b) in the case of the TD-STA method in the considered range of depths.For the TX8/RX32 data acquisition scheme the corresponding decrease in the contrast was approximately 11 dB (from −2 to −13 dB, Fig. 8a) in the case of the SD-STA method and approximately 5 dB (from −6 to −11 dB, Fig. 8b) in the case of the TD-STA, respectively.Also, slightly slower decrease of the CNR with growing depth was observed for images reconstructed using the conventional TD-STA as compared to the SD-STA.Specifically, the CNR decreased approximately by 9 dB (from 5 to −4 dB, Fig. 7a) in the case the SD-STA method and by approximately 3 dB (from 2 to −1 dB, Fig. 7b) in the case of the TD-STA method for the TX16/RX32 data acquisition scheme in the considered range of depths.The corresponding values for the TX8/RX32 data acquisition scheme were approximately 13 dB (from 5 to −8 dB, Fig. 7a) in the case of the SD-STA method and approximately 8 dB (from 0 to −8 dB, Fig. 7b) in the case of the TD-STA, respectively.
Comparing by visual assessment the B-mode crosssectional views (see Figs. 5 and 6 in Subsec.4.1) of the 3D numerical acoustic data simulated in Field II one can observe that the TD-STA method provided slightly better depth of visualization than the proposed SD-STA method especially for the cases of TX16/RX16 and TX8/RX8 data acquisition (see panels 4 and 5 in the Figs. 5 and 6).However, the SD-STA method provided better image resolution which can be visually assessed as image sharpness for the cysts over the entire range of depths considered.
Moreover, the proposed method allowed the time of the HRI synthesis to be reduced significantly due to the properties of the FFT algorithm (Brigham, 1988) utilized in the SD data processing.Specifically, the numerical complexity of the single LRI reconstruction was as large as O N x N y N t log(N x N y N t ) (Tasinkevych, 2017) for the case of the SD-STA method, where N x and N y are the number of samples in the synthesized image in the lateral plane and N tis the number of time samples recorded in a single RF acoustic echo.In the case of the conventional TD-STA method the corresponding value was O N x N y N z N 2 , where N z is the number of samples in the synthesized image in the axial direction, and N 2 is the number of RF signals used for the LRI synthesis.The corresponding reduction of the LRI reconstruction time was therefore (Tasinkevych, 2017): In the case of the TX16/RX32 data acquisition scheme (N = 32) and for the reconstructed images of 64 × 64 × 512 (N x × N y × N z ) samples where each RF had 1122 time samples (N t ) the corresponding value was: T SD−STA T TD−STA ∼ 1 20.This yielded about 80 times faster synthesis of the final HRI image comprised of 4 LRIs.
It is worth noting that to conduct a reliable comparative analysis of the results obtained using different methods (SD-STA and TD-STA) no data processing was applied to enhance the imaging quality, like dynamic apodization (Guenther, Walker, 2007;Mehdizadeh et al., 2012;Tasinkevych et al., 2013) or the time gain compensation taking into account the attenuation in the medium (Szabo, 2004), etc.This resulted among others in the relatively weak overall quality of the B-mode images obtained with the SD-STA and TD-STA methods, especially for the cases of TX16/RX16 and TX8/RX8 data acquisition schemes corresponding to the LRIs synthesis from the reduced matrices of RF echoes (see for instance Figs. 5 and 6 were the images obtained using the numerical data simulated for the set of anechoic cysts are shown).This examples clearly show the extent of final image quality deterioration due to acceleration of the acoustic data acquisition and processing.Therefore, further research is necessary in order to improve the quality of 3D ultrasound imaging with the use of the reduced RF matrices data acquisition schemes.Their development is especially important for the low-cost portable devices which should operate with a limited power consumption and will be the subject for future research.

Conclusions
In this work the SD-STA 3D image reconstruction method based on the non-overlapping TX/RX sub-aperture data acquisition combined with the spectrum domain data processing was proposed.For transmission of the interrogation signal and detection of the back-scattered echoes the limited number of signal channels (64 or 256) was used which determined the number of active elements utilized in the TX/RX events (8 × 8 and 16 × 16-element nonoverlapping TX/RX sub-apertures).Such a configuration of the ultrasound data acquisition schemes is consistent with the requirements of the low-cost portable devices usually operating within low power consumption limitation and using small number signal channels for transferring the acoustic data from the probe to the host.
The results obtained in this paper confirmed that a trade-off between the TX/RX sub-aperture size which determines the number of TX/RX events re-quired for the synthesis of the final image and its quality exists.Specifically, in this research the optimal subaperture size was 16 × 16-element one which required as much as 16 TX/RX events to reconstruct the HRI consisted of 4 LRIs, provided the full-size matrices of RF signals were used for their synthesis (TX16/RX32 data acquisition scheme).Applying the reduced RF matrix approach (TX16/RX16 data acquisition) the corresponding number of TX/RX events was reduced four times.However in the latter case the image quality deterioration was observed.
Moreover, the SD-STA method allowed a significant increase in the image reconstruction rate compared to the conventional TD-STA method.This acceleration was obtained at the cost of a slight deterioration in the image quality assessed by the contrast and CNR parameters which usually can be considered as tolerable in the case of the low-cost portable devices.
Further research will focus on improving the imaging quality.For this purpose the SD-STA method will be generalized for the case of overlapping sub-apertures.This however will inevitably lead to the frame rate decrease which is a crucial parameter in the realtime 3D ultrasonography.Therefore, optimization of the number and size of the TX/RX sub-apertures as well as their stride will be necessary to achieve the best possible quality maintaining the acceptable frame rate at the same time.Another possible way to improve the imaging quality is to generalize the SD-STA method for the case of the coherent compounding of ultrasound waves transmitted at different angles.In this case the numerical complexity of the SD data processing will increase significantly because one will have to deal with 3D interpolation of the spatial-temporal spectrum instead of the 1D interpolation which is required in the SD-STA method proposed in this paper.

Fig. 3 .
Fig. 3. Flow diagram showing the successive steps of the acoustic data processing in the considered SD-STA and TD-STA imaging methods.

Fig. 4 .
Fig. 4. Block diagram of the experimental set-up.Right panels show the real view of the scattering plastic sphere used in experimental measurements.

Fig. 5 .
Fig. 5. B-mode images of the numerical data simulated in Field II for the set of anechoic cysts.The cross-sectional view (plane Y = 0) of the 3D reconstruction obtained using the FD-STA method for different data acquisition schemes is demonstrated.images are plotted in logarithmic scale over 30 dB dynamic range.

Fig. 6 .Fig. 7 .Fig. 8 .
Fig. 6.B-mode images of the numerical data simulated in Field II for the set of anechoic cysts.The cross-sectional view (plane Y = 0) of the 3D reconstruction obtained using the time TD-STA method for different data acquisition schemes is demonstrated.The images are plotted in logarithmic scale over 30 dB dynamic range.

Fig. 9 .
Fig. 9. 3D images of the set of 5 scattering spheres reconstructed using proposed SD-STA method for different TX/RX data acquisition schemes.The images are visualized using sliceomatic visualization package for MATLAB in logarithmic scale.Isosurfaces corresponding to the −30 dB level are shown.

Fig. 10 .
Fig. 10.3D images of the experimental data acquired from the custom design scattering phantom reconstructed using the proposed SD-STA method for different TX/RX data acquisition schemes.The images are visualized using sliceomatic visualization package for MATLAB in logarithmic scale.Isosurfaces corresponding to the −30 dB level are shown.