Super-Resolution Ultrasound Imaging Using the Erythrocytes—Part II: Velocity Images

Super-resolution ultrasound imaging using the erythrocytes (SURE) has recently been introduced. The method uses erythrocytes as targets instead of fragile microbubbles (MBs). The abundance of erythrocyte scatterers makes it possible to acquire SURE data in just a few seconds compared with several minutes in ultrasound localization microscopy (ULM) using MBs. A high number of scatterers can reduce the acquisition time; however, the tracking of uncorrelated and high-density scatterers is quite challenging. This article hypothesizes that it is possible to detect and track erythrocytes as targets to obtain vascular flow images. A SURE tracking pipeline is used with modules for beamforming, recursive synthetic aperture (SA) imaging, motion estimation, echo canceling, peak detection, and recursive nearest-neighbor (NN) tracker. The SURE tracking pipeline is capable of distinguishing the flow direction and separating tubes of a simulated Field II phantom with 125–25-<inline-formula> <tex-math notation="LaTeX">$\mu \text { m}$ </tex-math></inline-formula> wall-to-wall tube distances, as well as a 3-D printed hydrogel micr-flow phantom with 100–60-<inline-formula> <tex-math notation="LaTeX">$\mu \text { m}$ </tex-math></inline-formula> wall-to-wall channel distances. The comparison of an in vivo SURE scan of a Sprague-Dawley rat kidney with ULM and micro-computed tomography (CT) scans with voxel sizes of 26.5 and <inline-formula> <tex-math notation="LaTeX">$5~\mu \text { m}$ </tex-math></inline-formula> demonstrated consistent findings. A microvascular structure composed of 16 vessels exhibited similarities across all imaging modalities. The flow direction and velocity profiles in the SURE scan were found to be concordant with those from ULM.


Highlights
• The paper hypothesizes that erythrocytes can be tracked in SURE to obtain vascular flow images and create velocity maps and profiles using a recursive nearest-neighbor tracker in recursive synthetic aperture ultrasound imaging.
• A Field II simulated phantom with 125-25-µm wall-to-wall tube distances, a 3-D printed hydrogel phantom with 100-60-µm wall-to-wall channel distances, and a Sprague-Dawley rat kidney are used for validation setup.
• The SURE tracking pipeline can distinguish the flow direction and separate tubes in the simulated and 3-D printed phantoms.The SURE velocity map has the same structure and flow direction as ULM and 26.5-and 5-µm microcomputed tomography (CT) scans.
into a variety of diseases, including cancer, diabetes, and other vascular diseases.
Several MB tracking methods are employed to reveal tracks and velocities including nearest-neighbor (NN) tracker [2], [17], bipartite graph-based method [18], modified Markov chain Monte Carlo data association method [8], Kalman tracker [10], [19], hierarchical Kalman tracker [4], and multifeature tracking [9].The MBs are injected into the bloodstream, and their density should be sparse to ensure separation for detecting and tracking individual bubbles to create a vascular flow image from the accumulated MB tracks.The sparsity of MBs requires an acquisition time between 30 s and 10 min to have enough MBs pass through the entire circulation to render a robust image.Furthermore, during a long acquisition time, precise motion correction and co-registration have to be made.This is difficult in a clinical setting because of both voluntary and involuntary movements.Additionally, MBs are fragile and may burst when exposed to a large acoustic pressure.As a result, the mechanical index (MI) of the sequence should be kept between 0.05 and 0.2, well below the current FDA limit of 1.9.
Contrast-free approaches for super-resolution imaging have been tried [20], [21], [22], [23], [24], [25].Bar-Zion et al. [24] used Doppler slicing for constructing microvessel images without contrast agents.They validated their method with two parallel blood vessel simulations and showed in vivo mouse brain images.It took approximately three minutes for a contrast-free super-resolved vascular image to be acquired with Doppler slicing, which is similar to a long acquisition of ULM with MBs. Park et al. [23] used deep learning to reconstruct microvessel images by localizing the positions of red blood cells.They used bifurcation phantoms and the great saphenous vein as an easy and nondense vascular structure without motion and ground-truth validation.You et al. [25], used a deep learning network based on power Doppler to reconstruct microvessel images.Their network was dependent upon ULM data and based on MBs for training.The dependency on their network for reconstructing microvessel images did not allow tracking scatterers to demonstrate vascular flow images and create velocity maps and profiles.A major disadvantage of all these three contrast-free super-resolution imaging methods was their inability to track scatterers and show vascular flow images, velocity images, and profiles since tracking is an integral part of super-resolution ultrasound imaging.These problems can be overcome without the use of contrast agents by using super resolution ultrasound imaging using the erythrocytes (SURE) [20], [21], [22].With approximately five million erythrocytes per mm 3 blood transporting oxygen throughout the body, there is a huge number of targets to detect and track in SURE.Therefore, SURE can be performed in just a few seconds of acquisition by accumulating detected erythrocyte peaks.SURE density imaging is presented in the accompanying paper [20].These papers demonstrate the use of erythrocytes for microvascular imaging.The main purpose of this article is to develop tracking methods to enable the quantification of flow velocities in the microvasculature.This entails compensating for tissue motion, echo canceling, and tracking to estimate the velocities.Tracking uncorrelated and high-density scatterers can be challenging and recursive imaging has to be employed to make it possible.The correlation between the detected positions in SURE is quite short, as isolated targets are not used.Fast imaging with low sidelobes is therefore needed to follow the flow.Low sidelobes are obtained by employing an interleaved synthetic aperture (SA) sequence with 12 virtual sources [26], and high frame rates are attained by using recursive imaging [27].Tracks can be made using an NN tracker to quantify flow metrics.The first work on tracking erythrocytes with an NN tracker was published in [17].Recursive SA imaging was introduced to track erythrocytes as high-density scatterers in a Field II simulated phantom, a 3-D printed phantoms, and a rat kidney [28], [29].This article expands on this work with a comprehensive version of the SURE tracking pipeline incorporating recursive SA imaging and recursive NN tracker.Field II simulated phantom is expanded with different tube distances in the same and opposite direction, a 3-D hydrogel printed phantom of varying channel distances, and an in vivo rat kidney, which is used to compare SURE with ULM using MBs, as well as a 26.5-µm voxel size micro-CT scan and a 5-µm voxel size micro-CT scan of the same kidney.
In Section II, the SURE tracking pipeline and evaluation methods are presented.In Section III, experimental results of a simulated phantom, a 3-D printed hydrogel microflow phantom, and a Sprague-Dawley rat kidney are presented.Then, the SURE image is compared with corresponding ULM and micro-CT images.Finally, a discussion and conclusion of this article is given.

II. SURE TRACKING PIPELINE AND EVALUATION METHODS
The SURE tracking pipeline is expanded from the pipeline introduced in the accompanying paper [20] describing the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.A. SURE Tracking Pipeline 1) Data Acquisition and GPU Beamforming: This article used three types of data, a simulated phantom using Field II Pro [30], [31], [32], a 3-D printed hydrogel microflow phantom, and in vivo data from a Sprague-Dawley rat kidney.The ultrasound data were acquired using a Verasonics Vantage 256 research scanner (Verasonics, Inc., Kirkland, WA, USA).All measurements were performed using a 168 elements 10-MHz GE L8-18iD (GE HealthCare, Chicago, IL, USA) lin-ear array probe.The pulsewidth modulation transmits stage in the Verasonics scanner yielded a center frequency in the tissue of 9.29 MHz and a wavelength of 166 µm, when estimated from the mean frequency of the received tissue signal's power density spectrum [20].For the SURE scan sequence, a pulse repetition frequency of 5 kHz with a high-resolution image (HRI) frame rate of 416.7 Hz was used.Full information and discussion about the SURE pipeline and scan sequence can be found in [20].Processing was carried out in MATLAB 2021b (MathWorks, Natick, MA, USA) with the processing modules explicitly developed for the SURE tracking pipeline.An interleaved SA sequence with 12 virtual sources with positive and negative transmissions was used [26].The stored data were beamformed using the graphics processing unit (GPU) code presented in [33].
2) Recursive Imaging: This article uses recursive SA processing [27], [28], [29].There are roughly five million erythrocyte scatterers per mm 3 in SURE compared to the sparse distribution of MBs in ULM.Therefore, erythrocytes have a high density in consecutive frames, making them impossible to track individually.It is, thus, the peaks in the speckle pattern, which was tracked.These peaks fairly quickly decorrelate due to the velocity gradients in the flow, and fast imaging is therefore needed.Here, the recursive imaging method can increase the frame rate 12 times due to the 12 virtual source transmissions to improve the tracking of the highly dense scatterers between consecutive frames.Schematics of conventional SA imaging and recursive SA imaging are shown in Fig. 2(a) and (b).In conventional SA imaging, an HRI is created by summing the 12 beamformed low-resolution images (LRIs), where each LRI is created after emitting with different transducer elements.When recursive imaging is used, a new HRI is created by summing the 12 beamformed LRIs after each new emission.HRIs are calculated recursively using the last 12 emissions, with the newest emission replacing the oldest.This gives 12 independent streams of frames for tracking the detected peaks after motion compensation and echo canceling.With this method, the frame rate increases from 416.7 to 5000 Hz, and there are 12 times more frames available for the detection and tracking of erythrocytes.
4) Erythrocytes Detection: The erythrocyte peaks are found from the remaining signals after echo canceling [20].Detection of erythrocyte locations is based on peaks above a specific threshold in decibels from the maximum of the image, and then, polynomial interpolation around the peak position determines the subsample location of each peak.To create the intensity map of SURE, this process is repeated for all frames in the sequence, resulting in a list of detected positions in all frames.Then, by having the list of detected erythrocyte positions, a velocity map can be created by tracking their positions.

5) Erythrocytes Tracking:
Tracking is an integral part of super-resolution ultrasound imaging, among them NN tracker is widely used to determine the closest targets in the next frame to the current targets [2], [17], [28], [29].Here, a recursive NN tracker is used to track high-density scatterers in SURE imaging.As mentioned in Section II-A2, recursive imaging updates the oldest beamformed LRI with a new one for creating an HRI.A schematic of the recursive NN tracker in recursive SA imaging is shown in Fig. 2(c).As a result of the 24 LRIs, recursive SA imaging produces 12 HRIs rather than two HRIs when conventional SA imaging is used.Here, rather than using the NN tracker for each HRI, 12 times NN trackers are used with 12 frame intervals in recursive imaging.The NN trackers are applied to HRIs with an interval of 12 as the 12 independent streams of data to insert more trajectories into velocity images.After this, all trajectories of these 12 NN trackers are summed up to create a velocity map.

B. Evaluation Method and Experimental Setup
This section presents the setup for the Field II simulated phantom, the 3-D printed hydrogel microflow phantom measurements, and the in vivo rat kidney measurements.
1) Simulation: A phantom with eight tube pairs was simulated using Field II [30], [31], using the Pro version [32], with 175 scatterers per resolution cell.This phantom consisted of eight pairs of tubes with axial center-to-center distances from top to bottom varying from 150 µm (0.90λ ) to 50 µm (0.30λ ) as listed in Table I.The tubes' radii were 12.5 µm, so the axial wall-to-wall distances of the tubes from top to bottom varied from 125 µm (0.75λ ) to 25 µm (0.15λ ).The last six wall-to-wall distance pairs of tubes and the last four center-center distance pairs of tubes were less than half of a wavelength (λ = 166 µm).Each pair of tubes had a parabolic flow profile with maximum peak velocities of ±10 mm/s.The first, third, fifth, and seventh pairs of tubes had opposite horizontal flow directions.However, the second, fourth, sixth, and eighth pairs had the flow in the  same horizontal flow directions.Using in vivo measurements obtained from a rat kidney [20], [28], a realistic tissue motion was applied to the simulated phantom.A frame rate of 416.7 Hz was obtained for 12 emissions, but after using recursive imaging, it increased to 5000 Hz.A total of 1250 frames were generated for three seconds of data.After recursive imaging, there were 15 000 frames.
2) Phantom Measurement: A 3-D printed hydrogel microflow phantom with horizontal microchannels was used for the experiment.A more detailed description of the fabrication process can be found in [36] and [37].Fig. 3 in the accompanying paper [20] shows the 3-D design of the channel structures and fiducial markers of the 3-D printed hydrogel phantom, as well as photographs of the sub-merged 3-D printed hydrogel phantom fixated between guide posts at its corners and of the measurement setup.As listed in Table II, the four channels' wall-to-wall separation from top to bottom varied from 100 µm (0.60λ ) to 60 µm (0.36λ ).The last two channels wall-to-wall distance was less than half of the wavelength (λ = 166 µm).Because of printer limitations, it was not possible to make wall-to-wall tubes closer than 60 µm.However, due to printing resolution and the hydrogel structure of the phantom, the channel sizes might be slightly different than designed [36], [37].The phantom channels were used with an infusion of highly concentrated MBs of SonoVue (Bracco, Milan, Italy) with a dilution of 1:1 and a flow rate of 0.2 µL/s.The channels' radius was 100 µm.The mean velocity of flow inside the 3-D printed phantom was calculated by dividing the flow rate of the syringe pump by the channel cross-sectional area, and assuming that the maximum velocity at the channel center to be twice the mean velocity as predicted for laminar Hagen-Poiseuille flow in a channel of a circular cross section.Therefore, the flow rate of 0.2 µL/s was equivalent to a maximum velocity of 12.7 mm/s.
3) Rat Kidney In Vivo Measurement: The in vivo experiments were conducted on healthy male Sprague-Dawley rats using the SURE imaging sequence [20].Here, 20 s of data were acquired with the SURE sequence on the rat kidney as described in Section II-A1.Also, 32 s of ULM with MBs was acquired with a frame rate of 139 Hz.An amplitude modulation SA sequence with 12 virtual sources was used for ULM, in which each virtual source emitted three times with all elements, even elements, and odd elements for full and half amplitude emissions.Both SURE and ULM imaging were performed with the same GE L8-18iD transducer at the same position of the rat kidney.Also, after SURE and ULM imaging, the rat kidney was prepared for the infusion of the intravascular micro-CT contrast agent Microfil (MV122, Flow Tech Inc., Carver, MA, USA) as a vascular reference image for comparison and validation.Ex vivo micro-CT imaging was performed for 11 h with a voxel size of 26.5 µm.In another scan, a smaller region of the rat kidney was scanned for 20 h with a voxel size of 5 µm.The goal was to identify the same structure in three images and determine to what extent there may be similarities and differences between the two methods of super-resolution ultrasound imaging [20].

III. RESULTS
Section III present the results of different experiments conducted with the SURE tracking pipeline.

A. Simulated Phantom Using Field II
The SURE resulting images from 3 s of simulated data are shown in Fig. 3.The B-mode image of one frame is shown on the top left of Fig. 3. Accumulation of detected peaks and trajectories map of the recursive NN tracker are shown as intensity and velocity maps on the top of Fig. 3.By using recursive imaging, 1250 frames within 3 s with a frame rate of 416.7 Hz was increased to 15 000 frames with a frame rate of 5000 Hz.The recursive NN tracker is applied to the erythrocytes' peak positions to create the SURE velocity map.The color wheel at the top right corner of the velocity map indicates the direction of the flow.Movement to the right or left is illustrated by the yellow or blue in the velocity map.The velocity map shows the first, third, fifth, and seventh pairs of tubes have flow in opposite directions; however, the second, fourth, sixth, and eighth pairs of tubes have flow in the same directions.The velocity profile of the tubes is shown in the lower Fig. 3.The black line is the estimated velocity profile, the green lines show the true velocities and the red lines show the velocities based on the theoretical model of a 33% underestimation of peak flow velocity in 2-D super-resolution ultrasound imaging according to [38] and [39].The estimated center-to-center distances of tubes, distance errors, estimated peak velocities, and velocities underestimation are listed in Table I.The estimated distance error for all pairs of tubes was in the range of 0-4 µm.According to Table I, the mean estimated velocities of the pairs of tubes in the opposite flow direction was 6 ± 0.8 mm/s (mean ± standard deviation), and for the tubes with the same flow direction was 6.4 ± 0.4 mm/s (mean ± standard deviation).The mean of the underestimation of peak velocities was 39% ± 8% for pairs of tubes in the opposite direction and 35% ± 3% for pairs of tubes in the same direction.

B. 3-D Printed Hydrogel Microflow Phantom
The SURE tracking pipeline was applied to three seconds of data acquired from a 3-D printed hydrogel microflow phantom.By using recursive imaging, 1250 frames within three seconds with a frame rate of 416.7 Hz was increased to 15 000 frames with a frame rate of 5000 Hz.Because there was no sample motion during scanning the printed phantom, the motion estimation stage was skipped in the SURE tracking pipeline in Fig. 1  2-D super-resolution ultrasound imaging [38], [39].As listed in Table II, the predicted center-to-center distances of channels are 305, 298, 289, 274, and 265 µm giving distance errors of 5, 8, 9, 4, and 5 µm.The wall-to-wall distances of the last two channels are 70 and 60 µm which is less than half a wavelength (λ = 166 µm).As a result of the velocity map and profiles of the 3-D printed phantom, all five channels can be correctly separated into two opposite flow directions.The errors between the estimated and actual channel distances were 5, 8, 9, 4, and 5 µm.The mean and standard deviation of the estimated peak velocities of five phantom channels were 8.6 ± 0.5 mm/s, giving an underestimation of 32% ± 4% from the true velocity of 12.7 mm/s (0.2 µL/s).

C. Effect of Recursive Imaging in the Tracking of High-Density Scatterers
In this section, the efficiency of recursive SA imaging and the recursive NN tracker were compared with conventional SA imaging.An 1-s data of the 3-D printed phantom from Section III-B was used with recursive SA imaging and conventional SA imaging.In conventional SA imaging, the NN tracker was applied for each HRI, which was created with 12 LRIs.For recursive SA imaging and recursive NN tracker, as discussed in Sections II-A2 and II-A5, HRIs were calculated recursively using the last 12 emissions, with the newest emission replacing the oldest emission.The 12 times NN trackers were used for the 12 different streams of HRIs from recursive imaging.The velocity maps and profiles for recursive and conventional SA imaging are shown in Fig. 5. Recursive SA imaging significantly increased the number of peaks and trajectories.In 1 s, the recursive SA imaging had 325 × 10 3 peaks and 9 × 10 3 tracks, but the conventional SA imaging had 14 × 10 3 peaks and 5 × 10 3 tracks.In the right graph in Fig. 5, the red line is one line of the axial cross-section velocity profile (the right white line in the velocity maps) of one second of data.The black line is the mean of 50 cross-section velocity profiles across a 50-µm width (left white line in the velocity maps).The standard deviation of the 50 cross-section velocity profiles is shown with the purple shaded area.The standard deviation of the 50 cross-section velocity profiles (left white line in the velocity maps) is shown with the purple shaded area in the velocity profiles in Fig. 5.The green dashed curves are the theoretical profiles.

D. SURE, Micro-CT, and ULM Comparison
The SURE and ULM images shown in Fig. 6 were acquired in 3 and 32 s, respectively.The 26.5-µm micro-CT scan was Fig. 5. Comparison of recursive NN tracker in recursive imaging and NN tracker in conventional SA imaging of the 3-D printed phantom for only 1 s of data.On the top is the conventional SA tracking image and its velocity profile.On the bottom is the recursive SA tracking image and its velocity profile.The red line is one line of the cross-section velocity profile from velocity maps (the right white line).The black line is the mean of 50 cross-section velocity profiles for a 50-µm width (left white line) from the velocity maps.The green dashed curves are the theoretical profiles.The standard deviation of the 50 cross-section velocity profiles is shown with the purple shaded area.acquired in 11 h and the 5-µm micro-CT scan was acquired in 20 h.For both the ULM and SURE images, the NN tracker was employed to track the locations of the observed peaks.By using recursive imaging in SURE imaging, 1250 frames within 3 s with a frame rate of 416.7 Hz was increased to 15 000 frames with a frame rate of 5000 Hz.The recursive NN tracker is applied to the erythrocytes' peak positions to create the SURE velocity map.The corresponding micro-CT images were matched to the SURE and ULM images.The SURE, ULM, 26.5-µm micro-CT, and 5-µm micro-CT from the same Sprague-Dawley rat kidney in the same region are shown in Fig. 6.The SURE intensity map had the same vessel structure as the 26.5-and 5-µm micro-CT image.The 16 markers in the images indicate the branches of the vessel structure in SURE and in the micro-CT images.These markers were also shown in the ULM image.The markers in the micro-CT, SURE intensity map, SURE velocity map, and ULM velocity map can be clearly correlated.The vessels markers numbers 13 and 15 were not clearly visible in the 26.5-µmCT image and were very thin in the 5-µm micro-CT image, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 7. From zoomed area 1.On the top from the left are the micro-CT image, SURE velocity map, and ULM velocity map (zoomed area of Fig. 6).The color wheel indicates the direction of the flow.On the bottom are velocity profiles of SURE and ULM from the white line indicated in the velocity maps.
but can be clearly seen in the SURE image and also in the ULM.
On the top of Fig. 7, a zoomed area of the 5-µm micro-CT image, the SURE velocity map, and the ULM velocity map are shown.The flow direction is the same in both the SURE and ULM velocity maps.A velocity profile of the same cross section of both SURE and ULM velocity maps is shown at the bottom of Fig. 7.The velocity profiles were interpolated and smoothed for better visualization.The estimated peak velocities of the first, second, and third vessels of the SURE image are 2.1, 1.1, and 2.5 mm/s and 1.8, 1.5, and 2.4 mm/s for the ULM.The vessel center-to-center distances of the first and second vessels in SURE are 236, and 243 µm for ULM.The vessel center-to-center distances of the second and third vessels in SURE are 269, and 290 µm for ULM.
On the top of Fig. 8, another zoomed area of the 5-µm micro-CT image, the SURE velocity map, and the ULM velocity map are shown.In both the SURE and ULM velocity maps, the flow direction is the same, but the structure of the SURE velocity map is more similar to the micro-CT image structure.A velocity profile of the same cross section of both the SURE and ULM velocity maps is shown at the bottom of Fig. 8.The velocity profiles were interpolated and smoothed for better visualization.The estimated peak velocity of the vessel of the SURE image is 2.5 and 1.8 mm/s for ULM.The vessel's diameter in the SURE image is 72 and 77 µm for the ULM, and 75 µm for the 5 µm micro-CT image, which are all less than half of the wavelength.

E. SURE Reliability
The reliability of SURE imaging and tracking can be revealed by making the images for independent data of the same object.This is attained by dividing the long dataset into segments of 3 s and then making the SURE density and tracking images for different segments.This is shown in Fig. 9 with starting times of 0, 6, and 12 s to ensure independent images.A video has also been made with a starting time increasing in steps of 0.25 s yielding 78 frames when the starting time is increased from 0 to 19.5 s.The video can be found in the supplementary material.

IV. DISCUSSION
The accompanying paper [20] demonstrated how to make the SURE density images by accumulating the detected erythrocyte peaks.Therefore, with SURE imaging, it is possible to obtain a super-resolution image without contrast agents in just a few seconds and with the full FDA range for MI.Results presented in this article demonstrate that erythrocytes detected in SURE imaging can be tracked to create a velocity map and reveal the velocity within vessels, enabling the separation of adjacent vessels.
In the simulated phantom, it was possible to separate the tubes with wall-to-wall distances of 125 µm (0.75λ ), 75 µm (0.45λ ), 45 µm (0.27λ ), and 25 µm (0.15λ ).The SURE tracking pipeline separated the eight pairs of tubes shown in Fig. 3 with the flow in the same and opposite directions.The wall-to-wall distances of the second, third, and fourth pairs of tubes were less than half of the wavelength, where the mean estimated center-to-center distances error were 2.3 µm in the opposite flow direction and 1.3 µm for the same flow direction.The results of Fig. 3 and Table I show that the SURE tracking pipeline could separate and estimate the flow direction not only for the wall-to-wall distance of less than half of the wavelength but also for the tube distances of 45 µm (0.27λ ) in both the Fig. 8. From zoomed area 2. On the top from the left are the micro-CT image, SURE velocity map, and ULM velocity map (zoomed area of Fig. 6).The color wheel indicates the direction of the flow.On the bottom are velocity profiles of SURE and ULM from the white line indicated in the velocity maps.
same and opposite flow directions.For a wall-to-wall distance of 25 µm, which is equivalent to 0.15λ , only opposite flow directions could be distinguished.It's important to note that this distance is three times smaller than the diffraction limit, which is half of the wavelength.The estimated distance error for all pairs of tubes was in the range of 0-4 µm, which indicated how accurate the SURE pipeline was in estimating the center-to-center distances of tubes with flow in the same and opposite directions.The mean of the underestimation of peak velocities was 37% ± 6% for all pairs of tubes in the same and opposite flow directions, which corresponds to the theoretical model of a 33% velocity underestimation [38], [39].According to the theoretical model introduced in [38] and [39], the peak velocity of a 3-D parabolic velocity profile is underestimated by up to 33% in 2-D imaging due to the limited spatial resolution in the elevation direction.In this case, the tube diameter in the simulated phantom was 25 µm, which was less than the 770 µm of full-width at half-maximum (FWHM y ) of the probe, i.e., the width of the vessel was much smaller than the FWHM in the elevation direction.There were some inaccuracies in the estimated velocities and tube distances for tubes with wall-to-wall distances less than 0.3λ , since the distances were too close compared with the third and fourth pairs with distances less than half of the wavelength (λ = 166 µm).
The 3-D printed hydrogel microflow phantom was designed and printed to investigate SURE imaging.To simulate high-density scatterers of erythrocytes in SURE, highly concentrated MBs were infused into the 3-D printed phantom.It was impossible to use Doppler fluid as blood-mimicking fluid due to stickiness within the channels of the 3-D printed phantom, leading to clogging.Therefore, the SonoVue was used to reuse the phantom multiple times.This does not mean that the SURE tracking pipeline is based on MBs; it only uses MBs as high-density scatterers in this particular experiment.It is not possible to print the channels closer than a wall-to-wall distance of 60 µm, which is still less than half of the wavelength (λ = 166 µm).Also, there is no specific criterion for measuring wall-to-wall distance from SURE intensity images.A reasonable way is to measure the center-to-center distance from the SURE intensity map.It is also possible that the distances between 3-D printed channels differ from those of the designed phantom.These errors may be because of shape deformation during the experiments, printing resolution, shrinking, or expansion of the phantom [36], [37].Using the SURE tracking pipeline with the recursive SA and recursive NN tracker provided accurate estimates of flow direction to the left, right, and down according to the velocity map shown in Fig. 4. It was demonstrated in Fig. 4 and Table II that the SURE tracking pipeline could separate the 3-D printed phantom channels with opposite flow directions.The wall-to-wall distances of the last two channels of the phantom were 70 and 60 µm, which is less than half of the wavelength (λ = 166 µm).The mean of the underestimation of peak velocities was 32% ± 4% for all pairs, which corresponds to theoretical model with 33% velocity underestimation [38], [39].The velocity maps for the channels in Fig. 5 show that recursive SA yielded a better velocity map and more trajectories with the recursive NN tracker compared to conventional SA after just one second of data.The velocity profile of one axial line cross section for recursive SA and recursive NN tracker was smoother than conventional SA.The mean for 50 cross-section velocity profiles of the recursive SA was smoother and more similar to the parabolic profile compared with conventional SA.Furthermore, the standard deviation of these 50 lines for recursive SA was significantly lower than that for conventional SA.This experiment utilizes data for only one second to demonstrate the impact of the recursive SA and recursive NN trackers on the velocity map and profile.Thus, the velocity profile is not expected to be similar to a laminar profile, as has already been demonstrated in Fig. 4 for three seconds of data.recursive imaging SA provides better detection and tracking of dense, non-separable targets, as well as a better velocity profile in less time than conventional SA.
The SURE image was compared with micro-CT images.The 26.5-µm micro-CT was scanned after a few days of extracting the kidney from the rat and embedding it into paraffin after injection of the micro-CT contrast agent.The micro-CT images may be changed in shape because the rat kidney is extracted from the rat and put into paraffine causing shrinkage.The slice of the 3-D micro-CT image was matched to the 2-D SURE and ULM images.It was possible to see some smaller vessels more clearly in the SURE and ULM images that did not appear in the 26.5-µm micro-CT image.For this reason, the top area of the rat kidney was scanned with a 5-µm voxel size a few months after the 26.5-µm scan.Therefore, the 5-µm micro-CT scan may further have changed in shape compared with the 26.5-µm scan.In Fig. 6, the markers indicated the vessel positions in the different images.Green markers indicate vessels in all images.Blue markers indicate vessels only present clearly in ultrasound images and 5-µm micro-CT image, but barely visible in the 26.5-µm micro-CT image.In Fig. 6, 16 points were marked to compare the SURE, the ULM, the 26.5 µm, and the 5-µm micro-CT scans.The markers number 13 and 15 could not be clearly seen in the 26.5-µmCT image.The reason for this is that the rat was euthanized after SURE and ULM scanning and there was no flow in the arteries.Injection of the micro-CT contrast agent causes the veins to increase more in size than the arteries, possibly due to the high elasticity of the arteries and the large compliance of the veins.By increasing the volume and pressure on the arterial side, the contrast will tend to accumulate on the venous side [40].Thus, a micro-CT with a higher resolution like the 5 µm was needed to visualize the arteries better.However, the SURE vessel structure density and thickness were closer to the micro-CT image than the ULM.The SURE velocity map shows flow directions approximately based on 16 markers in the ULM velocity map.Thus, it is possible to image renal microvessels in vivo in just a few seconds by detecting and tracking erythrocytes.The vascular structure is verified using micro-CT images obtained from an excised kidney during 11 and 20 h of scanning.
The ULM with MBs was acquired for 32 s in this experiment using a Verasonic scanner, but typically, these scans take 10 min to complete with a BK scanner to get a sufficient amount of data [4].In this experiment, it was impossible to acquire ULM data for more than 32 s using the Verasonic scanner.To obtain the same image location for both ULM and SURE image acquisitions, both were scanned using the GE L8-18iD probe and the Verasonic scanner for a fair comparison.
Comparison of SURE and ULM velocity profiles in Figs.7 and 8 showed that the SURE tracking pipeline can yield similar results as ULM.According to Fig. 7, the maximum peak velocities of those three vessels of SURE and ULM varied by 0.3, 0.4, and 0.1 mm/s, indicating how close their estimated maximum peak velocities were to each other.The SURE and ULM velocity maps of Fig. 7, which show the flow directions of the vessels, indicate that the SURE tracking pipeline was capable of correctly estimating the direction of the flow, which is shown by the color wheel.In addition, the center-to-center vessel distances of SURE and ULM varied by 7 and 21 µm, which indicated the distances of both SURE and ULM were quite close to each other.In another example, with a different zoom area of the rat kidney shown in Fig. 8, the 3-s SURE velocity map structure was similar to the 5-µm micro-CT image and the 32-s ULM.The vessel sizes of SURE and ULM were 72 and 77 µm with just 5-and 2-µm difference from the 75-µm vessel diameter in the 5-µm micro-CT image, where all vessels sizes were less than half of the wavelength (λ = 166 µm).Estimated maximum peak velocities for SURE and ULM differed by 0.7 mm/s.ULM velocity map looks poor in comparison to the literature and state-of-the-art results, but it should be noted that this is 32 s of ULM data, as discussed earlier due to the Verasonic limitations to obtain data over a sufficiently large time.Moreover, in the ULM velocity map, an NN tracker is used to have a fair comparison with the recursive NN tracker which is used in the SURE velocity map.It is important to note that ULM velocity maps can also be obtained based on using Kalman tracker [10], [19], hierarchical Kalman tracker [4], and multifeature tracker [9].This is beyond the scope of this article as it is intended to demonstrate an initial implementation of a SURE tracking pipeline using recursive NN tracks to demonstrate that velocity maps can be obtained within a few seconds using erythrocytes.

V. CONCLUSION
SURE imaging makes it possible to track erythrocytes to obtain vascular flow images and create velocity maps and profiles in a few seconds.The SURE tracking pipeline was used with modules for beamforming, recursive SA imaging, motion estimation, image alignment, SVD for echo canceling, peak detection, and recursive tracking with the NN tracker.The conducted experiments in this article included three different experiments of the simulated phantom using Field II simulation, a 3-D printed hydrogel microflow phantom with horizontal channels, and an in vivo experiment on a Sprague-Dawley rat kidney.The SURE tracking pipeline was able to distinguish the direction of flow and separate the tubes of the simulated Field II phantom with a wall-to-wall distance of 45 µm in the same and opposite directions.As well as this, wall-to-wall distance channels of 60 µm with opposite flow directions were separated for a 3-D printed hydrogel microflow phantom, where these distances were less than half of the wavelength.The SURE image of a rat kidney was compared to a ULM image, a 26.5-µm micro-CT image, and a 5-µm micro-CT image.The revealed vessel structures and flow directions were similar between images.In the comparison of the velocity profiles of SURE and ULM, it was discovered that the SURE tracking pipeline could yield similar results as ULM.Using the SURE tracking pipeline, it was possible to create a consistent velocity map within a few seconds of data acquisition.He is currently a Postdoctoral Researcher and works on image processing and analysis from computed tomography and magnetic resonance imaging.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 1 .
Fig. 1.SURE tracking pipeline stages.(a) Various processing blocks in the SURE tracking pipeline.(b) Images from each stage of the SURE tracking pipeline in the three different experiments are shown for the Field II simulated phantom, the 3-D printed hydrogel microflow phantom, and the Sprague-Dawley rat kidney.Image of scanner used by permission of Verasonic.

Fig. 2 .
Fig. 2. Schematic of conventional and recursive SA imaging, and recursive NN tracker.(a) Conventional SA imaging.An HRI is created by summing the 12 beamformed LRIs, where each LRI is created after emitting with different transducer elements.(b) Recursive SA imaging.A new HRI is created by summing the 12 beamformed LRIs after each new emission.HRIs are calculated recursively using the last 12 emissions, with the newest emission replacing the oldest.(c) Recursive NN tracker applied to HRIs of recursive SA imaging.The NN trackers are applied to HRIs with an interval of 12 as the 12 independent streams of data to insert more trajectories into velocity images.

Fig. 3 .
Fig.3.Results of the simulated phantom using Field II simulation.On the top from the left is the B-mode image, the SURE intensity map, and the SURE velocity map.On the bottom is the SURE velocity profile.The black line is the estimated velocity profile.The green and red dashed lines are the true peak velocity and the peak velocity predicted by the theory.
. The B-mode image of one frame is shown on the top left of Fig. 4. Accumulation of detected peaks and trajectories map of the recursive NN tracker are shown as intensity and velocity maps on the top of Fig. 4. The estimated velocity profile of the four pairs of channels with their estimated center-to-center distances is shown at the bottom of Fig. 4. The black line is the estimated velocity profile, the green dashed lines show the true velocities, and the red dashed lines show the velocities based on the theoretical model of 33% underestimation of the peak flow velocity in Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 4 .
Fig.4.Results of 3-D printed microflow phantom.On the top from the left are the B-mode image, the SURE intensity map, and the SURE velocity map.On the bottom is the SURE velocity profile.The black line is the velocity profile.The green and red dashed lines are the true peak velocity and the peak velocity predicted by the theory.

Fig. 6 .
Fig. 6.Comparison of SURE, ULM, and micro-CT images.On the top from the left is shown the SURE intensity map, the SURE velocity map, and the ULM velocity map.The color wheel indicates the direction of the flow.On the bottom is shown the corresponding 26.5-and 5-µm micro-CT images.Markers are vessel indicators for comparison of the structure.

Fig. 9 .
Fig. 9. Showing SURE reliability of the intensity map (left) and velocity map (right) when using data for 3 s with different start times.A video of all start times from 0 to 19 s with steps of 0.25 s can be found in the supplementary video .

Iman
Taghavi (Member, IEEE) received the bachelor's and master's degrees in electrical engineering from the University of Isfahan, Isfahan, Iran, in 2016, and the Ph.D. degree in super-resolution ultrasound imaging from the Technical University of Denmark, Lyngby, Denmark, in 2022.He has successfully implemented several digital systems and developed compressed sensing radars, achieving a significant reduction to 2% of the Nyquist criteria.His recent work has focused on the development of an advanced signal-processing pipeline for tracking microbubbles inside the bloodstream for super-resolution imaging.He is currently with the Center for Fast Ultrasound Imaging, Department of Health Technology, Technical University of Denmark, as a Postdoctoral Researcher.His research interests include signal processing and its applications in radars, array sensors, and biomedical imaging.Mikkel Schou (Member, IEEE) received the M.Sc.and Ph.D. degrees in 3-D ultrasound perfusion imaging from the Technical University of Denmark, Lyngby, Denmark, in 2017 and 2020, respectively.His research focused on row-column transducer arrays for perfusion and super-resolution imaging.He is currently employed at BK Medical and GE Healthcare, Herlev, Denmark, as an Ultrasound Specialist.Sebastian Kazmarek Praesius (Graduate Student Member, IEEE) received the M.Sc.degree in mathematical modeling and computation from the Technical University of Denmark (DTU), Lyngby, Denmark, in 2022, where he is currently pursuing the Ph.D. degree with the Center for Fast Ultrasound Imaging, working on real-time 2-D and 3-D beamforming, motion estimation, and super resolution.His research interests include graphics processing unit (GPU) acceleration, deep learning, and computer graphics.Lauge Naur Hansen (Graduate Student Member, IEEE) received the M.Sc.degree in quantitative biology and disease modeling from the Technical University of Denmark (DTU), Lyngby, Denmark, in 2022, where he is currently pursuing the Ph.D. degree with the Center for Fast Ultrasound Imaging, working on the validation of super-resolution ultrasound imaging techniques by quantitative comparisons with high-resolution computed tomography scans.Nathalie Sarup Panduro received the master's degree in medicine from the University of Copenhagen, Copenhagen, Denmark, in 2015.She is currently pursuing the Ph.D. degree in clinical and preclinical super-resolution ultrasound imaging with the Department of Diagnostic Radiology, Rigshospitalet, Copenhagen.Sofie Bech Andersen received the master's degree in medicine from the University of Copenhagen, Copenhagen, Denmark, in 2015.She is currently pursuing the Ph.D. degree in preclinical ultrasound super-resolution imaging in a collaboration between the Department of Biomedical Sciences, University of Copenhagen; the Department of Radiology, Rigshospitalet, Copenhagen; the Center for Fast Ultrasound Imaging, Technical University of Denmark, Lyngby, Denmark; and BK Medical, Herlev, Denmark.She was concluding a year as a Resident (introductory position) in radiology with Rigshospitalet.Her research interests include the development of diagnostic imaging tools for use in the clinic.Stinne Byrholdt Søgaard received the Ph.D. degree in preclinical ultrasound super-resolution imaging in a collaboration between the Department of Biomedical Sciences, University of Copenhagen; the Department of Radiology, Rigshospitalet, Copenhagen; and the Center for Fast Ultrasound Imaging, Technical University of Denmark, Lyngby, Denmark, in 2023.Carsten Gundlach received the Ph.D. degree in physics from the University of Copenhagen, Copenhagen, Denmark, in 2006.He is currently a Senior Researcher with the Department of Physics, Technical University of Denmark, Lyngby, Denmark.Hans Martin Kjer received the Ph.D. degree from the Technical University of Denmark, Lyngby, Denmark, in 2015.

TABLE I EXPERIMENTAL
RESULTS FOR THE SIMULATED PHANTOM WITH THE MAXIMUM PEAK VELOCITIES OF ±10 mm/s IN THE SAME AND OPPOSITE FLOW DIRECTIONS USING FIELD II

TABLE II EXPERIMENTAL
RESULTS FOR THE 3-D PRINTED HYDROGEL MICROFLOW PHANTOM WITH A FLOW RATE OF 0.2 µL/s EQUIVALENT TO 12.7 mm/s