Ultrasound imaging is an indispensable tool in medical diagnosis,1 primarily due to its cost-effectiveness and unique real-time capability. Contrast-enhanced ultrasound (CEUS) imaging aims to assess vascular flow and tissue perfusion. This requires the intravenous injection of gas-filled microbubbles (MBs), which are efficient point scatterers of ultrasound.2,3 Their diameter, typically around 2 to 3 μm, allows passage through the entire vascular bed. Their nonlinear response to ultrasound has been deployed with specially designed ultrasound transmit pulse sequences to enhance their signals while also suppressing tissue signals.1 However, after more than 20 years of research, there are very few CEUS applications in the clinic that are supported by health services around the world. The most important application is in the differential diagnosis of liver lesions,4–7 whereas other applications in the abdomen,8–10 cardiovascular,11–14 or using targeted and therapeutic MBs15–19 are rather research oriented or have limited clinical use. The main reason for the limited use of CEUS is the high interobserver and intraobserver variability compared with that what is achieved by other imaging modalities (eg, magnetic resonance imaging [MRI] or computed tomography [CT]). Although ultrasound imaging provides comparable resolution to MRI or CT, typically 0.3 to 1 mm for abdominal applications, a number of factors including the equipment type and settings, the patient variability, and contrast material20 affect its performance. All of the previously mentioned items have a significant role in the increased uncertainty of CEUS-derived measurements of blood volume and flow that are rather not quantitative.21,22 Despite their limited temporal resolution, MRI and CT are often preferred.
A recent advance, with significant potential in medical imaging, is the application of particle localization and tracking methodologies to CEUS, often inspired by light microscopy.23–26 Despite their small size compared with the ultrasound wavelength (typically around 100 times smaller), most MBs have high scattering cross-section (SCS) and provide strong ultrasound scatter that can be recorded as a single point spread function (PSF) in the image.27–29 Thus, their localization is possible, and early in vivo images showed at least a 5- and up to 10-fold resolution improvement at diagnostic frequencies (6.5 MHz).24,30,31 The work used a small number of single MBs that can be detected and then tracked. Algorithms for particle detection and localization can be implemented in ultrasound data.23,24,30,32–34 The focus is on localization of single MBs either isolated30,35,36 or with enough separation between echoes within a group of multiple MBs, and often overlapping MB echoes, due to MBs located closely, are rejected.24
In practice, the image plane includes a large dispersion of vessel diameters that range from a few micrometers and up to millimeters. The MB concentration will be linearly related to the volume of these different-sized vessel groups, and consequently, under conditions of constant volume flow, the number of MBs that cross any given microvessel will be a few orders of magnitude lower than the feeder vessel, which is millimeters in diameter. An infusion of contrast agent that ensures a small number of MBs in any given frame is likely to ensure single MB scattering events that are easy to detect, localize, and track. However, given the large number of microvessels and the need for enough MBs to cross the entire vasculature available in the image, it would be necessary to acquire long duration video loops in the order of tens of minutes or possibly hours depending on the application, size of organ, or region of interest. Such data sets render examination procedures impractical and are likely to provide significant motion artifacts that may be several magnitudes of orders larger than the aimed resolution. To alleviate this problem, the length of video data is required to reduce to a few minutes and thus the concentration of MBs must be increased. Tracking of a large number of MBs is also possible.24,25,31,32,35 The most efficient method is provided by high frame rate imaging that can provide scatter overlap for flowing MBs and deploy tracking using a simple image analysis method that follows the MB path by means of autocorrelation.23,26 This means that all single MBs will be tracked and the tracking efficiency is near 100%. To achieve high frame rate, single-emission protocols, for instance plane waves, were used. This provides a high acoustic pressure that drops significantly with depth and normally a high acoustic pressure near the transducer, which has 2 consequences: (1) the penetration depth is limited as there is no focus to offset the attenuation, and (2) the MB detection sensitivity is highly variable. This is due to the highly variable acoustic pressure with depth and the dependence of the MB SCS on acoustic pressure27,37 (ie, SCS increases with increasing acoustic pressure). In general, higher acoustic pressures provide larger MB densities as a larger number of MBs scatter above noise. In addition, if the acoustic pressure is high then there will be significant MB destruction. Thus, the assumption that the MB density across the image is uniform and representative of the vascular volume is difficult to uphold.
A near-uniform pressure field can be achieved by using a focused transmission. It is possible to approximately offset the attenuation using a focus at the bottom of the image. Although the attenuation is variable at different parts of the image, this approximation ensures a near-uniform distribution of MBs across the image that provides a good approximation of the concentration of red blood cells and hence blood volume. Currently, ultrasound imaging systems use nearly exclusive focused transmissions, and it is also possible to generate ultrasound images with single MBs using these systems. However, the frame rate will be low in the order of 10 to 20 Hz. This makes MB tracking challenging as MB echoes across subsequent frames usually show no overlap, complicating the reconstruction of the MB trajectories. Also note that the speed of red blood cells ranges from micrometers to centimeters per second.38,39 Thus, in the case of large MB density in the image, it is difficult to identify each MB in consecutive frames. Furthermore, a large MB density will provide several MB events that are merged and not separated. Although the focus on isolating single MBs24 ensures good localization results and thus optimal resolution outcome, the merged MB events are data that should be used to maximize data usage and reduce video loop duration. More recently, sparse recovery methods seem promising in providing similar quality information by using large MB densities and deploying prior knowledge of the PSF.40,41 Here we use lower MB densities that would enable their separation but would also create several merged MB scatter events in the image. It is assumed that these merged MB events are likely to be in image regions of larger vessels where the flow is pulsatile, hence the MB density is high. On the other hand, microvessels may provide mostly single MB events as the blood volume drops and the velocity profiles are nearly close to an average.
A model-based approach for tracking a large number of particles has been developed for tracking single fluorescent molecules in optical microscopy in images that include underlying clutter.42,43 The tracking method is based on (a) a detection and segmentation algorithm that localizes the particles retaining their image features into a particle probability image (PPI), and (b) correlating their image features to detect their tracks.42–44 If tissue motion is not a problem, there is no limit to the duration of the video loop that this algorithm can process, and the longer the video loop, the more tracks will be detected. Tracks are formed by linking detected particles so their number depends on the number of detected particles per frame, which is assumed to be fixed per frame. Wilson et al43 showed good performance in images with low signal-to-noise ratio (SNR) in microscopy applications; these images contain items similar to CEUS level of noise. In this work, image sequences of approximately 2000 frames (30 frames per second) were processed, and over 4000 tracks were created (SNR, ~2 dB) from 25,000 detected particles. These image conditions resemble CEUS images with a sparse MB distribution. Thus, this method can be applied to CEUS. The model-based approach allows the inclusion of prior knowledge on the evolution of MB signals45–50 that may lead to more efficient tracking. Emphasis is given to the detection and segmentation. This is because, unlike optical microscopy, ultrasound images have a highly variable PSF (ie, an individual MB will have a different appearance in different locations within the image), which is a significant challenge for the detection process. Methodologically, a ground truth is required for testing the algorithms, and this is provided by simulated environments. To date, super-resolution ultrasound work provides little information on localization accuracy and particularly on the role of detection and segmentation in this. Most experiments use vessel environments that do not provide this information,23,25,51,52 but rather demonstrate that the vessels can be resolved and that the MBs remain within the vascular space. The most common measure of such an accuracy is the full-width at half maximum (FWHM) of the diameter of each vessel24,31 or the size of a single MB PSF.30 Here, the detection and localization accuracy are assessed using synthetic data that simulates real data from animal or human studies. The accurate assessment of the number of detected particles by this approach, including spurious and missed detections, and the localization uncertainty are important for both methodology development and measurement of performance. As a second step, an in vitro vascular phantom validation is required. These methods do not replace the synthetic data validation above as no information on MB location, other than the fact that they are located inside a vessel, is provided. Instead, the aim is to assess image resolution. In addition, and by using real imaging equipment at settings that are relevant to clinical imaging, a convincing case is presented that MBs are detected and tracked as they flow inside microvessels with a resolution beyond the diffraction limit. A well-controlled experimental setup with a translucent phantom was used32 here.
However, such a setup has low variations of speed of sound compared with those generated by tissue. Thus, this idealized setup provides lower aberration in the near field of the array and the PSF of the entire setup is less variable compared with clinical imaging. This presents less challenging conditions to the operation of the detection, localization, and tracking methodologies. For these reasons, we chose to use an in vivo validation by means of optical coherence tomography (OCT). This ensures a live in vivo validation with exactly the same equipment, imaging parameters, and contrast agent, with real vessels and tissue aberration. The animal model used is the ovine ovarian corpus luteum (CL), a highly vascularized gland, which has previously been shown to be suited for perfusion imaging studies.53 Two further ex vivo methods for the microscopic examination of tissue, namely, optical microscopy and optical projection tomography (OPT), have been deployed to identify vascular features across the entire tissue that compare well with the ultrasound images. Finally, a proof of concept examination of a prostate patient is presented alongside the ex vivo histological evaluation for comparison. Imaging tumor vascularization is an active research field,21,54–58 including prostate tumors.59,60 The study of prostate cancer diagnosis includes all imaging modalities, and MRI seems as the most promising.61
MATERIALS AND METHODS
The ultrasound image analysis was developed using synthetic, in vitro, and in vivo animal data. The latter used the sheep ovary and had 2 purposes: first, to provide the basis for generating synthetic data with a good approximation to in vivo image data; second, to validate the application of the methodology in vivo for clinically realistic ultrasound imaging parameters.
Animal Contrast Ultrasound Imaging
All animal data were collected under UK Home Office license approval (Duncan PPL 60/4401), and all the methods have been previously described.53 For each experiment, a female adult sheep was anesthetized, and each ovary in turn was exposed by laparotomy. The ovary was secured in place while maintaining full blood supply and was covered with a coupling ultrasound gel to achieve optimal contact with the ultrasound probe.53 This ensured that its shape was not distorted due to breathing motion and during imaging. Off-plane movements were minimized by using a reticulated arm that fixed the probe in a position where the scan plane was parallel to the breathing motion.
The methodology developed here is intended for use with standard CEUS that is currently available in the clinic. The contrast administration protocols were thus designed to provide a high density of single MB scattering events in the image, aiming to reduce the duration (eg, a few minutes), similar to that of a clinical patient examination. An iU22 Philips ultrasound scanner (Philips Medical Systems, Bothell, WA) with an L9–3 linear array probe (3 MHz) was used throughout (λ = 514 μm). Two contrast administration protocols were used. A 1.2-mL bolus of SonoVue (Bracco, Geneva, Switzerland) contrast agent and an infusion of 1.2 mL of SonoVue in 20 mL of saline at 0.5 to 1.5 mL/min rates using a Vueject pump (Bracco SpA, Geneva, Switzerland). All video data were saved in the DICOM format for offline processing.
In Vitro Imaging
To assess the performance of the super resolution image analysis in vitro, a cellulose tube, extracted from a single use renal dialysis filter (Dialyzer GFE-09, Gambro, Germany), secured to a 28-gauge microfill needle (World Precision Instruments, Stevenage, United Kingdom) with Luer connection was mounted in a Perspex tank. The tank was filled with boiled and cooled water. This setup was based on previously described in vitro work.47 The cellulose tube has a nominal internal diameter of 200 μm. Subsequently, and upon inspection using bright field microscopy, it was found to have a maximum internal diameter of 300 μm. The same ultrasound scanner and transducer as described previously were used. The transducer was mounted above the tube. 1.2 mL SonoVue contrast agent was added to 50 mL of water and allowed to flow by gravity feed through the cellulose tube. The end of the tube was placed out with the Perspex tank so no MBs entered the surrounding water. B-mode and contrast mode video loops of MBs flowing through the tube were acquired, saved, and later processed offline.
Tissue Processing and Microscopy Imaging
To obtain information on the vascular structure of the scanned ovaries, 3 different tissue imaging methods were used: standard optical microscopy, OPT, and OCT. The OCT was undertaken in vivo after CEUS images were acquired. For microscopy and OPT, the ovary was removed after CEUS and processed as described later. As the only live validation method, OCT was expected to provide the most accurate data on vessel diameter and structure. This validation method serves a dual purpose: (a) it advances on an in vitro validation because the in vivo setup here provides the additional challenge of real tissue imaging aberration due to the variable speed of sound provided, and (b) it can be used as the best criterion standard for comparison with the other 2 optical methods that are performed ex vivo. Note that it is very difficult to provide MB location data using an optical validation method, so this task is undertaken by using the synthetic data.
Optical Coherence Tomography
Live images of the structure and information on the dynamics of the vascular bed were acquired using an OCT with a Telesto-II (wavelength 1300 nm, 5.5-μm axial resolution in air; Thorlabs, Lübeck, Germany).62 The OCT probe was secured directly over the top of the ovary and perpendicular to the ultrasound probe. This arrangement aimed at capturing the same image plane with both methods. This was a challenging experimental setting as the ultrasound probe had to be horizontal, imaging the top 1 to 2 mm of the tissue. Apart from the suboptimal gel coupling, the positioning of the ultrasound probe was not informed by the live OCT image, as the microvessels that appeared in the OCT are impossible to visualize in B-mode. The OCT data comprised a 3-dimensional (3D) data set of the light scattered from the tissue and blood at the surface of the ovary. Optical coherence tomography files were saved in ThorImageOCT software (Thorlabs, Lübeck, Germany) and exported to ImageJ for analysis.63,64
For ovaries undergoing histological slicing and staining, the process has been previously described.53 The removed ovaries were fixed, dehydrated, and embedded in paraffin, after which 5-μm-thick sections were sliced. For each 10th section, biotinylated lectin BS-1 (Sigma, United Kingdom), used to specifically bind to endothelial cells, was applied overnight at 4°C in a humidified chamber. Impact diaminobenzidine (Vector Labs, United Kingdom) was used to stain the endothelial cells. Sections were assessed by light microscopy (Olympus BH2) and a digital camera; images were saved at ×20 magnification. Vessels across the whole plane could be measured and counted. For comparison with ultrasound imaging, the slice taken from a location as close as possible to the imaging plane was chosen.
Optical Projection Tomography
Previous work formed the basis of the method here.65 Before removal of the ovary, 70-μL rhodamine-labeled Griffonia (Banderiaea) simplicifolia lectin 1 (GSL 1 lectin; Vector Labs, United Kingdom) in 7 mL of phosphate-buffered saline was flushed through the ovary followed by 7 mL of phosphate-buffered saline, using the main feeding artery of the ovary as the input source. The GSL 1 lectin was used to stain the main vessel endothelial cell walls. The ovary was then removed and stored in 4% paraformaldehyde. Before OPT, the ovary was immersed in 1 part benzyl alcohol (Sigma Aldrich, United Kingdom) and 2 parts benzyl benzoate solution (Acros Organics, United Kingdom) for 2 to 3 days to achieve tissue translucency. A Bioptonics 3001 OPT scanner (Bioptonics, Edinburgh, United Kingdom) was used to acquire tomographic scans generating 801 image files per scan at a rotation of 0.45 degree. Acquired images were 1024 × 1024 pixels with pixel size of 18.28 μm. The ovary tissue was autofluorescent in the green channel, and the GSL 1 lectin staining was visible on the red channel. Scan planes were combined to form a 3D reconstruction of the ovary using NRecon software (Skyscan, Kontich, Belgium). Interrogation of the 3D volume to select the same plane as that of the ultrasound image plane was undertaken in Bioptonics 3D Volviewer and in ImageJ (NIH.gov, United States).
Synthetic Data Formation
The first step was to simulate the MB movement within a vessel network. The simulation provided the spatial coordinates and identity of each particle at any time. This information was used as the ground truth for testing the subsequent localization. Such ground truth information is not available in in vitro or in vivo experiments. The synthetic vessel network consisted of a 1.1 × 2.2 cm grid of interconnecting vessels with diameters that ranged from 10 to 500 μm, while their length ranged from 25 μm to 2 mm. The particle flow within that network followed Poiseuille's law: by applying first a constant global pressure drop across the network and second the mass conservation law, the pressure field could be calculated. When the calculated elemental pressures were substituted in Poiseuille's law, the volumetric flow rate within each vessel was obtained. The flow simulation method is illustrated by a set of diagrams that are included in Boujelben et al.66 The injected dimensionless particles moved according to the calculated flow value.
Ultrasound Image Simulation
The next step was to assign an MB echo to each particle position. It is known that the PSF is variable across the ultrasound image. As MBs are nonlinear scatterers of sound, they will provide echoes that are highly variable and difficult to model. To our knowledge, there is no research or commercial ultrasound image simulator that provides a reasonable MB echo simulation. Thus, we avoid a full calibration of the PSF here and instead generate a distribution of MB echo characteristics from in vivo video loops drawn from the sheep ovary experiment. A common characteristic of the MB echo image in this experiment was the ellipsoid shape of variable orientation (Fig. 1A). Three-dimensional Gaussian Fitting (minor, major axis, and intensity information) was used to extract the topology profile of the MBs in the in vivo data set. The resulting distribution provides a range of sizes, intensities, and shapes for MB echoes. This process can also be seen as an approximate estimation of the PSF. Note that it is accepted that the axial resolution is at best half the pulse length,67 which in our case is 1 λ (514 μm), whereas the lateral resolution varies across the image depending on the beam width67,68 and can be up to a few λ using the delay-and-sum beamfomers,69 which are used here. The histogram of the MB echo size distribution from the sheep ovary is displayed in Figure 2. The first and third quantiles are 8 (0.14 mm2) and 20 (0.35 mm2) pixels, respectively, whereas the median is at 12 pixels (0.21 mm2). Single MB echoes are unlikely to have sizes above 12 pixels (0.21 mm2), and the minimum detected MB is 5 pixels (0.09 mm2, approximately slightly larger than 2 × 2 pixels or one half λ × one half λ). Note that once 2 particles approach closer than a few wavelengths, they are likely to become larger in size before appearing as a merged particle in the image.69
In addition, the noise was characterized in the background and added to the simulated system as white Gaussian noise. This is because Rayleigh noise becomes additive Fisher-Tippet noise after log-transformation,70 which is well approximated by additive white Gaussian noise.71 Two hundred synthetic frames with simulated MBs were created. Figure 1B shows the result of the above process compared with a frame from a video loop of sheep ovary (Fig. 1A). Figure 1, C and D are the magnified squares from Figure 1, A and B, respectively, for further detail.
Super-Resolution Image Analysis
An MB tracking algorithm was developed here to process the ultrasound video loops. The algorithm structure was based on a single particle tracking software previously used in optical microscopy.43 This software follows four main steps: (a) to detect particles, (b) to segment them (ie, identify their region), (c) to localize them (ie, find their center or exact location), and (d) to track them as they move from frame to frame. This process ensures that a highly accurate particle location is identified (super-resolution localization) and then its path is generated by joining particle locations at consecutive frames. The result can be displayed in particle density, path density, or velocity maps that contain different structural and dynamic information. Here, the particles are MBs that move only within the lumen of blood vessels. Thus, their localization is a location inside a vessel, and their path identification is equivalent to drawing the vessel structure. Because ultrasound imaging has different and highly variable aberration compared with optical microscopy, its PSF is highly variable. This results in particle appearance that varies across the image. Thus, the detection and segmentation part of the algorithm is central to the development here to suit the ultrasound MB/particle application.
The CEUS images were processed to remove artifacts, such as bright echoes that were identified in the preinjection images, and were cropped from the entire frame sequence using ImageJ. Further, as off-plane movements were kept to a minimum, the well-established rigid registration7 was used to remove in-plane movement due to respiration. This was done using B-mode image features at the end-expiratory phase as previously described.72
A semiautomatic detection process was designed to detect image particles with a range of sizes using an adaptive nonlocal means filter,44 which has shown to work well with super-resolution ultrasound microvessel imaging.73 The manual input of 3 parameters was required: the average MB echo intensity, the minimum and the maximum MB echo size. These are roughly estimated in an initial observation of the video loop. The MB intensity was preserved by using the PPI,43 a nonlocal mapping of the original grayscale image through the use of multiscale Haar-like features.44 These features were obtained after the convolution of 3 kernels, with variable formation and size that was determined manually according to the range of the particle sizes. The resulting Haar feature image is the linear combination of the maximum value of each spatial scale at each pixel. In the normalized Haar feature image, the pixels with higher values considered as statistically significant, which means they are likely to belong to an MB. A noise threshold is implemented in this image, which is the average noise level of the original frame sequence. This results in a binary image that is used to calculate the PPI. Subsequently, the particle probability in each pixel is the ratio of the number of pixels that have value equal to one in the binary image divided by the total number of pixels of the particle area.
Furthermore, a region threshold of 1/e of the normalized PPI is applied to initially estimate the target regions.43 This enables the discrimination of the foreground and the background and generates an initial segmented image. In parallel, the input image is convolved with a Gaussian kernel to create a Gaussian smoothed image, which is used to extract the local maxima of the image. Only the local maxima points that belong to a segmented region and have particle probability value above the threshold are preserved creating the final watershed seed points. These points determine the regions that the final segmentation process is based on. The final detection is refined by an automatic update in the last step of the detection process. Microbubble echo regions with a size below the minimum input size and input average intensity and above the maximum input size, which depended on the SNR of the frame sequence, were eliminated. This serves as a third noise classifier as all particles below the minimum size and average intensity effectively are considered as noise. This process also enabled the automatic MB detection update from frame to frame.
Segmentation is the key step to accurate localization. Marker-controlled watershed segmentation, as described, works efficiently in optical microscopy, where the PSF of the system is known and is generally symmetric and fairly constant in a frame sequence. In live-cell imaging, the algorithm implements segmentation by using a dilated gradient image of the Gaussian smoothed image. This image and the gradient image are the input values of the watershed transform. Subsequently, the watershed transforms segments of each particle. The change in particle area that is due to the dilation process does not affect the localization accuracy as the intensity-weighted center of mass method works well in a circular and constant PSF.
However, in CEUS, the PSF is variable not only between different frames but also within the same frame as the acoustic field and aberrations change across the image. As a result, MB echoes have nonregular shapes, which make the accurate segmentation of each region essential for efficient localization. Thus, the gradient image was replaced by an inverted Gaussian as the input variable of the watershed transform. The latter can better recover the MB echo area and avoid the reduction caused by the gradient image. This proved particularly useful for MBs that have low intensity and thus low SNR as well as overlapping MBs that are close to each other and were difficult to discriminate.
As mentioned previously, overlapping MBs are a significant number of detected events in the ultrasound imaging. By including these in the processing, data usage is maximized and it is possible to generate maps with short video loop duration. In super-resolution ultrasound, often a number of MBs situated close to each other, creating a large common echo, were mostly eliminated, and only clearly discriminated single (individual) MBs were included in the detection process.2,3,24,25,30–32,35,52 In optical microscopy, the accurate knowledge of the PSF permits the differentiation of particles that partially overlap in the image,43,74 and all the particles are possible to localize. This is difficult to implement in CEUS images due to the lack of PSF constancy.
However, it is possible to treat overlapping MB events as single. The rationale is that a single center has high likelihood to be located inside its vessel. The overlapping MBs may follow the same detection and segmentation process as a single event. Thus, the maximum MB size is set to a large value to process regions of both single and overlapping MBs. The algorithm can determine the size and the neighborhood of the detected MBs by adjusting 2 parameters, the Gaussian smoothing and the local maxima width. Fewer detections may occur if either the Gaussian smoothing or the width are increased. Because in an overlapping event there may be low confidence as to whether the local maxima truly represent single MBs, the choice of larger values in the Gaussian width and the local maxima width enables a single detection. Finally, overlapping events follow the same process of the birth and death of single MBs within the algorithm, permitting the split and merge of the events,42 which may occur in consecutive frames.
Microbubble Localization and Tracking
The final segmented regions (MBs) were used for the localization of each MB. The intensity-weighted center of mass43 that is commonly used in ultrasound field24,30,33 was also used here. This is because the final segmented regions preserve all the original information to deploy this method. Alternative methods include the use the FWHM,25 a deconvolution of the PSF,23,26,31 or the local maxima to calculate the center of the PSF. These methods do not deploy the intensity of the segmented regions or rely heavily on the assumption that the PSF is constant and were not the first choice here. Microbubble linking in consecutive frames (ie, tracking) was completed using the nearest neighbor method.42
Following the linking process, all the tracks were superimposed in one figure, creating this way the final density map. The pixel value in this map was the number of tracks that passed through this pixel. Further, this algorithm provided information about the speed of the MBs, and an MB velocity map was generated, which represented the mean velocity of the tracks that passed through each pixel. The pixel size can be at the original size or at a subdivision of the original size reflecting the localization uncertainty. Through synthetic data experiment, the level of the subdivision was determined by creating a subpixel similar to the root mean square error (RMSE) of the localization uncertainty.
The detection efficiency and localization accuracy of the algorithm was tested by means of five statistical measurements on the synthetic data as commonly used in optical microscopy.75 These are as follows:
- Root Mean Square Error. The RMSE evaluates the localization accuracy of the software using paired particle events. Each ground truth event is paired with the result of the algorithm that is in its vicinity. The distance between the 2 provides a measure of the localization deviation distance. For the number of pairs in the image sequence, this is given by the following formula:
where et is the localization deviation distance for the pair t. Thus, the lower the RMSE value, the better the localization.
- Missed Events. Any ground truth event that does not pair with a detection event is counted as a missed event.
- Spurious Events. Any detected event that has not been attributed to a ground truth event is counted as spurious events.
- Minimum Distance. This is the minimum localization deviation among paired events.
- Maximum Distance. This is the maximum localization deviation among paired events.
In Vitro Data
The algorithm performance was tested by comparing the resulting vessel diameters on the final density maps with the diameter of the tube.
In Vivo Data
The accuracy and efficiency of the methodology was evaluated by comparing features on the final density maps with those identifiable on optical microscopy, OPT, or OCT. Cross-section area measurements from the CL and the follicles as well as vessel diameters from arteries, arterioles, and follicle wall vessel diameters were used.
Proof of Principle Prostate Patient Data
The method was also tested in vivo for the ability to resolve the prostate microvasculature in a patient referred for radical prostatectomy at the AMC University Hospital (Amsterdam, the Netherlands) because of biopsy-proven prostate cancer. An intravenous peripheral injection of a 2.4-mL SonoVue bolus was performed. The bolus passage through the prostate was recorded for over 2 minutes by transrectal ultrasound imaging with an iU22 ultrasound scanner (Philips Healthcare, Bothell, WA) equipped with a transrectal probe C10-3v. For the acquisition, a power modulation pulse scheme at 3.5 MHz was adopted to achieve contrast-specific imaging, and the mechanical index was set equal to 0.06 to avoid bubble destruction. The frame rate was equal to 10 Hz. The probe was held at a fixed position, and no breathing movement was observed. After radical prostatectomy, the resected prostate was first fixated and then cut in slices of 4-mm thickness, which were then marked by a pathologist to delineate prostate cancer based on microscopic analysis of cellular differentiation.
Development of the Super-Resolution Algorithm on Synthetic Data
A real-frame sequence from a sheep ovary was used to generate a sequence of synthetic ultrasound frames. Typical frames are shown in Figure 1, A and B, respectively.
Figure 1, D to F illustrate the result of the new segmentation process. The original synthetic image (Fig. 1D) was segmented using the watershed function. The result from implementing 2 different methods, the gradient image, as used in optical microscopy (Fig. 1E), and the inverted Gaussian (Fig. 1F) are shown. As a first observation, the size of the PSF of the segmented MBs were smaller (number of segmented pixels displayed next to each MB in Figure 1D) when the gradient image was used (Fig. 1E) compared with that of the inverted Gaussian (Fig. 1F), which had very good agreement with the original data (Fig. 1D). This slight difference between these 2 sets of measurements (Fig. 1, D and F) is based on the existence of the noise in Figure 1D and the manual calculation of the size of each MB there, whereas in Figure 1F, the same calculation took place after the process of noise removal by the algorithm.
Table 1 summarizes the difference between the 2 segmentation approaches. Better size and shape recovery (Fig. 1F) of the original MBs (Fig. 1D) led to the significant improvement in localization accuracy. The RMSE was 3 times lower using the inverted Gaussian (25.8 μm) compared with the gradient image (78.6 μm). Second, there was a significant decrease of the missed events using the inverted Gaussian method to 0.07% (Table 1). The high occurrence of missed events using the Gradient image is explained by the removal of scatterers by the algorithm, which was significantly reduced in size using the gradient image (arrow in Fig. 1, E and F) as they have size below the MB echo input size. As a result in Figure 1E, 2 of 5 detections are recorded, whereas in Figure 1F, 5 of 5.
Overlapping events where MBs were merged in appearance in the synthetic data were assessed. An example is shown in Figure 1G where 2 MBs are 59 μm apart. If the algorithm is used to deploy the detection of local maxima to identify each of the 2 MBs (Fig. 1G), they are (Fig. 1H) 396 μm apart, and 172 μm and 224 μm away from the tube, respectively. Increasing the typical particle size that reduces the segmented regions, the 2 merged MBs were treated like a single event (Fig. 1I). The distance of this event from the center of the tube was 14 μm. The processed synthetic video loop contains 10,780 single MBs in 200 frames, which resulted in 2442 overlapping MBs (mostly due to double particle overlap) and 5896 well-separated single MBs (8338 in total). The overlapping events were effectively 2 or more MBs that were in a distance below the three fourth of the wavelength (385 μm) from each other. These were treated as single events. The number of both single and overlapping events is in a good correspondence with the software result that detected 8117 events.
Assessment of Super-Resolution Methodology in Real Imaging
Performance Assessment In Vitro
The results for the in vitro feasibility test are shown in Figure 3. The video loop comprised 300 frames and was collected using the same ultrasound settings as the in vivo work. A B-mode frame and the corresponding contrast mode frame are displayed in Figure 3, A and B, respectively. The methodology using the inverted Gaussian image was used to process the in vitro frame sequence. This was limited to an ROI that was away from artifacts (yellow box of Fig. 3B). Figure 3C presents an example processing with detections within the red box of Figure 3B. The events that result from segmenting overlapping events into 2 or more single ones (green circles) and those that result from treating overlapping events as single (red crosses) are presented. It can be seen that in treating overlapping events as single leads to localizations inside the tube while segmenting these into single events leads to 2 localizations that are 600 μm apart (tube maximum internal diameter, 300 μm). Thus, segmenting overlapping events is often erroneous and leads to an overestimation of vessel size.69 Further, the density maps at a 3 × 3 subdivision of the original pixel size were studied using the 3 available detection processes. Figure 3D shows the density map that resulted from detections of only single and well-separated events, whereas Figure 3E used single events that stem after the segmentation of overlapping events, and Figure 3F detects overlapping events as single along with the well-separated single ones. The FWHM at the same 3 example positions (1, 2 and 3) was measured for each of the above processes and are given in Table 2. This shows that the diameter estimation is best depicted by the map of Figure 3F, where overlapping events are treated as single. Note that the lack of data in Figure 3D leads to an underestimation of the diameter and the lack of accurate PSF knowledge leads to an overestimation of the diameter in Figure 3F. Finally, Figure 3H shows the corresponding density map for Figure 3E using the original pixel size.
The vessel diameter is depicted accurately when overlapping MBs are treated as single and confirms our hypothesis on the overlapping MB location. These results are in agreement with the synthetic data. By treating overlapping events as a single event, all localizations are found inside the vessel, and the vessel width is depicted accurately. On the other hand, our knowledge on the PSF is not adequate to inform the splitting of overlapping events accurately (similar to Fig. 1H), and therefore they can be located outside the vessel resulting in an overestimation of the vessel width (Fig. 3E). In addition, when only well-isolated single MBs are used, the number of detections drops by 40%, compared with when overlapping events are accounted for (Table 2), which leads to an underestimation of the vessel diameter.
Algorithm Performance Assessment In Vivo
The extracted ultrasound video loops comprised 500 to 1500 (processed) frames with frame rate of 12 to 13 Hz. All videos were acquired in the wash-in or the wash-out period of a bolus injection apart from the videos used for OCT comparison, which were generated by using an infusion of sparse MBs. The mean video duration for bolus injections was 121 ± 19 seconds and for infusions 92 ± 28 seconds.
Figure 4 illustrates an example of the algorithm performance in an in vivo setting depicting the density maps using the original pixel size (132 μm) or the 3 × 3 subdivision of the original size (44 μm) as described previously. Figure 4A shows the tiled microscopy image, lectin stained to show endothelial cells, of a slice of a typical mature CL of the sheep ovary (day 9). The circular pattern is due to the near-spherical shape of the CL. This shape is the result of the rapid angiogenic growth of the CL in the first 9 to 10 days of the oestrous cycle of the ewe. Its vascular architecture here shows that the largest arterial feed surrounds the tissue and branches into smaller arterioles. The size of this artery is between 200 and 300 μm. The arrows in Figure 4A show that the artery is no longer in view in this slice, which is due to continuation and branching off plane. Branching of this artery into arterioles is evident across the entire length of these arteries. The core of the CL is mostly filled with microvessels. Not shown in Figure 4A is the supporting ovarian vessels (artery and vein) that lie under the location of the CL. Figure 4, B and C are the B-mode and peak contrast image in CEUS mode, respectively. The MB movement was tracked through the CL, 27 seconds after the peak contrast image and when the population of MBs was sparse enough for the algorithm to operate. The video loop consists of 583 frames, and its duration was 48.6 seconds. In Figure 4D, the processing includes the inverted Gaussian, and the parameters are adjusted for single MB detection only. The number of detections were 10,906, the average number of detections per frame was 19 (ranging from 1 to 42), and the number of tracks was 1135. Note that, in the absence of ground truth, the number of missed or spurious events cannot be quantified. Figure 4E displays the path density map that uses the inverted Gaussian in the segmentation process and is set to detect both single and overlapping MBs. This resulted in 34,325 detected events, which is on average 59 events per frame (ranging from 7 to 138), and 2951 tracks. The display has a pixel size 44 μm. Figure 4F shows the corresponding mean velocity map. Figure 4G shows the density map for single and overlapping MBs just as in Figure 4E, but uses the gradient image in the segmentation process instead of the inverted Gaussian. Finally, although Figure 4, D to G used the 3 × 3 subdivision process, the respective density map of Figure 4E with the original pixel size is shown in Figure 4H.
There is good resemblance between Figure 4A and Figure 4E, whereas Figures 4D and 3G depict a number of sparse tracks that bear little resemblance to the vascular architecture of the sheep ovary (Fig. 4A). Arrows in Figure 4E show a very similar location of the termination of the feeder arteries (arrows in Fig. 4A). Further, these larger vessels have a large number of paths (in yellow) compared with the inner CL, which agrees with their large comparative size shown in the histology. In addition, a number of smaller vessels are shown to branch inwards from these feeder vessels as shown in the histology slice. The rest of the inner CL area has a small number of paths that may be attributed to low microvessel flow. Under the boxed area of Figure 4E is the ovarian artery and vein, not shown in Figure 4A. These display the largest density of paths in Figure 4E and largest blood velocity (Fig. 4F). These features are not apparent in Figure 4D and Figure 4G that use only single MB processing and gradient image approach in the segmentation, respectively. In addition, using the original pixel size for localization in density maps (Fig. 4H) cannot depict different types of vessels. The feeder vessels that surround the CL (200 to 300 μm) and other microvessels inside the CL do not appear different to the ovarian artery and vein below the CL that are a couple of mm wide.
Density maps were obtained and processed similarly from 6 contrast video loops, which were typical of the range of data sets that were acquired, and the information is displayed in Table 3. From these, 1 was in wash-in and wash-out of the bolus, limiting the process only within the CL, 3 were in the wash-out and 2 were captured during an infusion of sparse MBs. The duration of each video loop varies from 39 to 129 seconds.
Table 3 shows that detecting both single and overlapping MBs using the new methodology with the inverted Gaussian as the input for the segmentation process maximized the number of detected events. As there is no ground truth in these data, a manual observation confirmed that the proportion of missed (undetected MBs) and spurious events (wrong detections) is low and did not affect the final density map. Indeed, within the CL, there is clear vascular network pattern. This table shows that the overlapping and single MBs using the inverted Gaussian (column 7) was always the maximum number of detected events in each data set and thus optimal detection. When the algorithm was adjusted to count the single events (column 6), it was shown that the overlapping events were between 32% and 68% of the total number of detected MB events. Columns 4 and 5 show the particle underestimation effect of using the gradient image as input for the watershed function instead. When both single and overlapping events were used in the detection process, the number of events accounted for were between 17% and 57% less than the ones that were counted using the inverted Gaussian. This was due to the shrinking of these particles' areas resulting from the gradient image calculation. The shrinking caused a lot of MB events to diminish to sizes that were not more than 5 pixels, which is the minimum input size threshold, and thus were eliminated. In addition, the size discrimination between single and overlapping events became more difficult to implement. This is because the algorithm size threshold was implemented after the segmentation process, and the shrinking due to the gradient image resulted in a lot of overlapping MB events to be misclassified as single. The resulting column 5 is likely to include mostly overlapping events, as many single ones were removed by the minimum input size threshold. In addition, when the particle size parameter was adjusted for single events, then the number further decreased, but a large number of overlapping events remained included (column 4). If it is assumed that the total number of true MB events is approximately that of column 7 and that the single events are reasonably approximated in column 6, it is then possible to produce an estimated number of single MB events accounted for in column 4, which is the difference between the column 5 number and the number of overlapping events (difference of columns 7 to column 6). This varies between 0% and 63% of the total number calculated. Negative numbers in brackets signifies that no single events are accounted for, and a number of overlapping events were also missed. These results confirm the synthetic data behavior observed in Figure 2, D and F.
Validation of In Vivo Ultrasound Size Measurements
Density maps from 10 different ovaries from 8 different sheep were compared with either histology (4 ovaries), OPT (6 ovaries), or OCT (2 ovaries) imaging. From these, 8 were bolus injections and 2 were captured during an infusion. Optical coherence tomography provided a live image of the vessels close to the surface of the ovary in situ, and the comparison with the respective density map is displayed in Figure 5. The main vessels are clearly seen in the density map (Fig. 5A) and the OCT image (Fig. 5B). Table 4 provides a summary of the measurements made. For the images shown in Figure 5, the vessels were measured to have diameters of 0.9 ± 0.01 mm from OCT compared with 1.19 ± 0.01 mm on the density map, and the narrow vessel diameter was measured to be 228 ± 23 μm from OCT compared with 236 ± 15 μm. The mean difference in measured sizes between the 2 was approximately 10%, which is not significant. Thus, OCT provided the best agreement between density maps, showing that the density map (Fig. 5A) can provide quantitative and accurate information of vascular structures, including vessel diameter and bending.
Compared with the OCT, the OPT provided information from larger areas of the tissue and in 3D. Compared with standard optical microscopy, OPT provided a 3D reconstruction of the whole ovary. It was then possible to choose a subvolume at a region of choice, which could also have a similar thickness to the ultrasound scan plane, which was approximately 2 mm, for comparison. This is not possible with standard 2-dimensional (2D) microscopy where the slices have a thickness of approximately 5 μm. Their orientation is roughly estimated before slicing, and therefore the comparison with the ultrasound image is difficult to make. Thus, in the example of Figure 4, it is rather fortuitous that Figure 4A appears very similar to Figure 4E. An example of a density map and the corresponding OPT slice are shown in Figure 6. The structure of the ovary from the ultrasound image can be seen in each modality: (a) peak contrast in CEUS mode, (b) density map in super-resolution processing, and (c) OPT. There are structures in the ovary identifiable in both the density map and the OPT slice such as the follicles (F1, F2), the CL, and larger vessels (V). The 2 follicles are well defined by the density map, and the large CL has a dense vascular network. On the smaller scale in the density map, there is a track detected in the first follicle, which is assumed to be part of the outer surface of the follicle. This is also seen in the OPT image (arrowed). This vessel measures 141 ± 54 μm diameter on the density map compared with 83 ± 7 μm on the OPT.
The main comparison performed between density maps and histology slices was on the measurement of the area of the CL. For 3 density maps and histology pairs, the mean CL area measured 48% smaller on histology than that measured on the density map. For all measurements, the mean difference in sizes measured on OPT is 36% smaller than that measured on the density map. Of the measurements made on the density maps, the smallest measurable regions that could be compared with similar regions on OPT were small vessels and follicle walls. A vessel diameter measurement on OPT at 83 ± 7 μm is compared with 141 ± 54 μm on the density map. Follicle wall thickness was measured to be 173 ± 35 μm on OPT compared with 241 ± 10 μm on the density map. It is evident that histology and OPT provide a significant underestimation of all sizes and area measurements compared with the density maps.
Finding comparable vessels in the optical images that were smaller than 100 μm and could be compared with the same vessels in the density maps was challenging. This is because all the different criterion standard techniques did not provide an abundance of such measurements. The OPT staining was inadequate to delineate arterioles, whereas the OCT resolution and sensitivity are not adequate to show the smallest vessels. As mentioned previously, the 2D microscopy with lectin staining provides such vessels but cannot be matched with the ultrasound maps. A number of vessels below 100 μm were observed, however, in the ultrasound maps and are displayed in Table 5 but were not verified. However, the large number of tracks strongly suggests the existence of vessels. The narrowest vessel measured on the density maps was 55 ± 10 μm, and there are several clear vessel paths between 50 and 100 μm.
Example Map of the Prostate
The MB tracking algorithm was applied to a human prostate with cancer. Figure 7, A and B are the B-mode and a contrast image after the peak of the MB density in CEUS mode, respectively. The acquired data set was typical and of low SNR. The video duration was 136 seconds from which the 115 seconds was the processed time that provided a density and velocity map (Fig. 7). The detection parameter combination was optimized to detect both single and overlapping MBs. A total of 1149 frames were processed, where 612,439 events were detected and created 47,536 tracks. The examination time had to be kept as short as possible; hence, an MB bolus injection was used. As a result, there is high MB density per frame, where the average number of detections per frame was 533.
The marked areas by a rectangle shape, an oval shape, and an arrow in Figure 7, C and D show the tumor areas that correspond to the histology. In the histology (Fig. 7E), cancerous areas are displayed with red color, based on the microscopic analysis of cell differentiation (Gleason's pattern) of whole-mount sections according to Montironi et al.76 Before the histopathological analysis, the prostate was fixated in formalin for 24 hours and sectioned in 4-mm slices. Note that the histological slices and ultrasound imaging planes are not parallel, and one imaging plane can intercept multiple slices. The velocity map provides very good correspondence to the histology particularly slices 5, 6, and 7, whereas the density map provides also a reasonable agreement with these slices.
Super-resolution images under realistic patient imaging conditions were achieved, demonstrating the feasibility of clinical 2D ultrasound super-resolution imaging using a standard CEUS mode. The gain in resolution is at least 5-fold, as vessels under 100 μm were detected at transmit frequency of 3 MHz (λ = 514 μm), and the system resolution here is approximately λ (half the pulse length). The smallest verified vessel width was 60 μm (Table 5), and the unverified detection of small arterioles (55 μm) presented nearly an order of magnitude resolution gain. The synthetic data investigation shows that the MB localization uncertainty can achieve 26 μm accuracy. The use of synthetic data enabled the development of the method into an ultrasound one as the errors induced by several parts of the processing were possible to assess and minimize. This was done by investigating detection efficiency, segmentation accuracy, and subsequently MB localization accuracy. The in vivo results are comparable with the literature in terms of resolution improvement. Experiments in thinned skull of rats in a fixed location provided λ/6 resolution at 20-MHz transmit frequency and using ultrafast scanning.31 Elsewhere, λ/4 resolution was achieved using higher transmit frequencies for identification of tumors.33 A clinical scanner has been used in the initial demonstration of super-resolution imaging in vivo, and under 20 μm resolution was achieved, which is over 5 times the improvement to the system resolution.24 This was performed in an optimal setting to minimize aberration, as only a thin slice of tissue was scanned, a flattened mouse ear, and the depth was also limited to 1 cm, and the tissue and probe were static. Our results that resolved structures with at least λ/8.5 accuracy were performed under conditions that are closer to clinical 2D CEUS using standard CEUS mode, low frame rate, radiology-relevant image depth to investigate a volume of tissue. Thus, compared with the above studies, significantly increased aberration was present in the in vivo data here. This results in additional PSF shape distortion with potentially negative consequences in MB detection and segmentation. Given the approximately 2-mm thickness of the 2D ultrasound plane, it seems surprising that vessels with tens of micrometers in diameter can be visualized. This is explained as the tracking algorithm enables the super-resolved MB localization across the third dimension of the slice thickness. The resulting density and velocity maps are in fact a projection of the ultrasound-exposed 3D volume into a 2D image plane. Thus, provided that scattering events are possible to distinguish, vessel structure is preserved projected and so is blood velocity. This is because different vessels may lie at different angles (ie, from vertical to parallel) in relation to the scan plane. This may result in inaccurate depiction of velocity as the angle is not known. However, in principle, 2D super-resolution imaging in vivo is not hampered by the width of the ultrasound scan slice, and several microvessels are possible to depict in density maps. Although it can be argued that in the future the velocity accuracy may be improved using 3D CEUS, the 2D prostate image here (Fig. 7D) strongly suggests that the velocity accuracy is not a problem as the high blood velocity were detected only around the tumor and correlated well with its area. This may be attributed to the high density of tumor neovessels that also have irregular pattern. This ensures that a large number of vessels are parallel to the scan plane. All this suggests that, in a clinical study, different types of velocity maps need to be tested to identify those that represent closest tumor dynamics in 2D (eg, maximum velocity that may be hypothesized to represent parallel vessel velocities). In addition, tissue motion artifacts do not appear in the super-resolution literature and seem to be well compensated here. Note that the animal experimental setup ensured that nonrigid motion is avoided and that the rigid motion due to breathing is kept in plane with the 2D ultrasound image plane. Thus, the well-established rigid registration provided good compensation. Further, the prostate did not require tissue motion compensation as breathing motion does not affect the position of the tissue.
The short video loop time (approximately 2 minutes) here, which provided adequate data for processing, strongly suggests that a clinically relevant examination time is feasible. However, such reduction of data results in additional challenges for the processing compared with other studies. The short video loop time required the use of a large number of MBs per frame. Table 3 shows the number of the detected events (single-overlap column) and the time for each processing, giving an overview of the average number of detections per frame. Comparing with literature, our algorithm detects, for the case of prostate cancer, 533 events on average per frame in 1149 frames, whereas the corresponding values, for example, in Errico et al,31 are approximately 13 events per frame, for a 75,000-frame data set. As mentioned previously, millimeter-sized vessels have orders of magnitude more blood volume, and thus MB concentration, compared with microvessels. At sparse MB infusion concentrations, this, in theory, results in very few single MBs in the microvascular bed, whereas a lot more and several overlapping MBs will appear in the larger vessels. Indeed, it was found here that more than 32% of detected events are attributed to overlapping MBs, which implies that overall these account for more than 50% of the MBs in the image, as shown in the CL study (Table 4). Previous investigations tend to avoid overlapping MB events.24,31,32 In this case, large data sets are required to depict the entire vascular space under investigation,24,25,31,52 which implies that clinical examination times would be significantly increased. The advantage of excluding overlapping events is that no assumptions are needed to include these events, and the localization accuracy is optimized. However, using only the single events the visualization of large vessels in their entirety may not be depicted accurately. The comparison between Figure 4D (single only) and Figure 4E (overlapping MBs included) showed that, within the approximately 49 seconds video loop time, the inclusion of only single MBs in the processing provided maps that do not include larger arterioles and feeder vessels. In other words, the exclusion of overlapping events provides a systematic error in mapping the vascular bed. The inclusion of overlapping MB echoes in the processing provided a more accurate depiction of the vascular structure with better representation of MB path proportion in different-sized vessels. It is suggested that the path density correlates well with volume flow, whereas the exclusion of overlapping MBs results in a blood volume estimation error for larger vessels. This further strongly suggests that the assumption that most overlapping events are likely to be located in the larger vessels is correct.
The challenge of including large MB numbers in the imaging may be best addressed using high frame rate imaging, which can provide MB scatter overlap and deploy tracking using the autocorrelation method.31 As mentioned in the introduction, such frame rates require a plane wave transmission that provides high attenuation and limited penetration. Further in CEUS, this limitation further increases the variability of the MB detection efficiency across the image due to increased S/N variability and MB destruction. Thus, the detected MB density and path density do not correlate well to red blood cell density and volume flow, respectively, which severely limits quantitative super-resolution maps of vascular dynamics. Here it is demonstrated that the less variable field of the focused transmission at low nondestructive acoustic pressures ensures reasonably uniform MB detection, with high penetration depth up to at least 6 cm (Fig. 7A). In addition, and as mentioned previously in Errico et al,31 an average of 13 events per frame were detected, whereas the method presented here processed, in the case of the prostate (Fig. 7), enables the handling of a larger number of detections as more than 700 events were detected in some frames. The tracking, used here, used a combination of the nearest neighbor approach and knowledge on the MB intensity, suggesting that it is not significantly inferior to the autocorrelation method. The tracking is in effect a sparse recovery method for the location of an MB path that has very few MB detections. In addition, there is no evidence in the literature that suggests that a high frame rate improves the statistics of the processing. Given that a bolus injection requires a minimum of a couple of minutes for the first pass to complete in most organs, it is suggested that it is this MB transit time, along with the dimensions of the vascular bed, that determine the MB number that is adequate to map the entire vascular structure. Here, it is suggested that mapping capillaries may be beyond the capability of image-based methods that aim at high resolution. Thus, vessels of tens of micrometers in diameter are the realistic targets of super-resolution ultrasound in clinical radiology, and most of these may be crossed several times in by the MB concentrations used here, thus providing adequate signal for processing. The inclusion of overlapping MBs is thus necessary. An approach that deals with the cumulative signal from all MBs in each pixel rather than particle events and, thus, uses the localization of spatially nonisolated MBs77 may be argued to include single and overlapping events in the processing. However, it is not evident from this processing what are the MB event areas, and thus it is much more difficult to assess detection and localization accuracies.
The super-resolution density maps presented here are the result of a robust detection and segmentation processes. All of the MBs were used in the tracking process including both single MBs of low intensity and overlapping MB events of large size and intensity. The signal enhancement through PPI and Haar-like features and the noise removal discriminating the background from the foreground make this algorithm capable of detecting even the weakest scatter. Christensen et al24 based the detection of single MBs on the cross correlation of each region in the reference frame and the subsequent frames. The groups that use ultrafast imaging detect individual events based on the correlation of each MB with the corresponding one in the next consecutive frame, due the high frame rate.23,26,31,51 These approaches are not effectively different to the one proposed here. Our approach enables the automatic detection and identification of MB events. Other methods apply the detection in B-mode frames after a similar approach to our background and foreground discrimination,32,35 thus dealing with multiple MB numbers. However, in CEUS, the shape of each scatterer in the image is nonregular and needs to be preserved. As a result, an optimal segmentation, using appropriately the watershed function, may be proposed for accurate localization of each event as well as ensure that nearly all MBs are detected.
Initial images of prostate cancer have been presented here from one patient. Although not conclusive, this initial result seems promising. Both density and velocity maps show good correlation with the histological evaluation. The velocity map suggests that tumor areas have redundant anastomotic vessels due to neovascularization. Extensive research aims at detecting and grading cancer by imaging technology to replace the use of invasive systematic biopsies (standard procedure) and reduce the risk of overtreatment and undertreatment. As aggressive (high grade) prostate cancer is correlated with angiogenesis and increased microvascular density,78,79 the proposed method may represent an asset for improved prostate cancer diagnostics and monitoring.
In future work, immunohistology by, for example, CD31 staining should be used to establish a ground truth of the microvascular architecture and improve over the adopted indirect comparison with the histology. A full study is required to assess this information that is otherwise difficult to compare with histological evaluation. The different microscopic techniques provided a criterion standard. However, these are limited, and this has become apparent through the improved near-microscopic resolution performance of ultrasound super-resolution images. Both microscopy histology and OPT images provided significantly reduced size measurements compared with the ultrasound density maps (Table 5; 41% mean difference), which confirms several other studies. It is known that optical imaging is undertaken ex vivo and after further tissue processing. This processing results in tissue shrinking.80 Further, fixation and histological preparation distort the tissues, and some of the variables required for this are difficult to standardize.81,82
The OCT performed in vivo provided the best comparison of vessel sizes with the density map with measurements within 10%. This confirmed that the measurements in the density maps are very accurate and demonstrated that ex vivo criterion standard comparison may be less appropriate for super-resolution ultrasound development. In addition, this OCT validation is superior to that using an in vitro setup with a capillary or capillary network that has been used previously.23,30,32 The in vitro setup provides lower PSF variability across the image compared with in vivo tissue imaging. This is because significant changes in the speed of sound across tissue affect the aberration and augment PSF variability compared with that of a translucent in vitro setup. This may be the reason that in vitro setups have been used previously to validate velocity estimation,32 as the diameter estimation is not a robust validation for the real imaging in vivo. In this context, the vessel diameter and thus system resolution are better validated in vivo. The 2 different-sized vessels (Fig. 5), which were embedded in tissue, had different vessel diameter and curvature, which provided a convincing validation of the vessel diameter accuracy of the ultrasound density maps at multiple positions in the same image. However, the most significant challenge with the use of the OCT was the matching of its plane with that of the ultrasound image. First, the ultrasound image is of low resolution and thus impossible to compare in situ. Second, the OCT is limited in depth (1.5 mm), which imposed an unusual position for the ultrasound transducer and restricted the imaging plane to very close to the surface of the ovary (and perpendicular to the plane in view for the rest of the imaging). However, it was shown that live techniques are more likely to be of use in the development of similar high-resolution ultrasound-based imaging methods.
In conclusion, a new super-resolution tool, which can be used with current clinical 2D CEUS, was presented. The feasibility was demonstrated in vivo for clinical radiology relevant image depths, with tissue motion and for short examination times under 2 minutes. The potential lies in the identification of regions with abnormal vasculature and particularly malignant tumors.
We also acknowledge the support of Sebastian Schaefer, Thorlabs (Lübeck, Germany) for the loan of and help with Optical Coherence Tomography; Harris Morrison of the IGMM, University of Edinburgh, for invaluable support, advice, and help with Optical Projection Tomography; Joan Docherty for facilitating and animal management for the ultrasound work; Nadia Silva for histology and experiment management; and Linda Nicol for the advice, ovary washing, and fixing.
1. Hoskins PR, Martin K, Τhrush A. Diagnostic Ultrasound Physics and Equipment
. 2nd ed. Cambridge, United Kingdom: Cambridge University Press; 2010.
2. Sboros V. Response of contrast agents to ultrasound. Adv Drug Deliv Rev
3. Sboros V, Tang MX. The assessment of microvascular flow and tissue perfusion using ultrasound imaging. Proc Inst Mech Eng H
4. Claudon M, Dietrich CF, Choi BI, et al.; World Federation for Ultrasound in Medicine; European Federation of Societies for Ultrasound. Guidelines and good clinical practice recommendations for Contrast Enhanced Ultrasound (CEUS) in the liver - update 2012: a WFUMB-EFSUMB initiative in cooperation with representatives of AFSUMB, AIUM, ASUM, FLAUS and ICUS. Ultrasound Med Biol
5. Dai Y, Chen MH, Yin SS, et al. Focal liver lesions: can SonoVue-enhanced ultrasound be used to differentiate malignant from benign lesions? Invest Radiol
6. Bartolotta TV, Taibbi A, Midiri M, et al. Indeterminate focal liver lesions incidentally discovered at gray-scale US: role of contrast-enhanced sonography. Invest Radiol
7. Ta CN, Eghtedari M, Mattrey RF, et al. 2-tier in-plane motion correction and out-of-plane motion filtering for contrast-enhanced ultrasound. Invest Radiol
8. Chung YE, Kim KW. Contrast-enhanced ultrasonography: advance and current status in abdominal imaging. Ultrasonography
9. Huang DY, Yusuf GT, Daneshi M, et al. Contrast-enhanced ultrasound (CEUS) in abdominal intervention. Abdom Radiol (NY)
10. Brabrand K, de Lange C, Emblem KE, et al. Contrast-enhanced ultrasound identifies reduced overall and regional renal perfusion during global hypoxia in piglets. Invest Radiol
11. Ten Kate GL, van den Oord SC, Sijbrands EJ, et al. Current status and future developments of contrast-enhanced ultrasound of carotid atherosclerosis. J Vasc Surg
12. Appis AW, Tracy MJ, Feinstein SB. Update on the safety and efficacy of commercial ultrasound contrast agents in cardiac applications. Echo Res Pract
13. Schinkel AF, Kaspar M, Staub D. Contrast-enhanced ultrasound: clinical applications in patients with atherosclerosis. Int J Cardiovasc Imaging
14. Schneider M, Anantharam B, Arditi M, et al. BR38, a new ultrasound blood pool agent. Invest Radiol
15. Gao S, Zhang Y, Wu J, et al. Improvements in cerebral blood flow and recanalization rates with transcranial diagnostic ultrasound and intravenous microbubbles after acute cerebral emboli. Invest Radiol
16. Wang S, Wang CY, Unnikrishnan S, et al. Optical verification of microbubble response to acoustic radiation force in large vessels with in vivo results. Invest Radiol
17. Porter TR, Xie F, Lof J, et al. The thrombolytic effect of diagnostic ultrasound–induced microbubble cavitation in acute carotid thromboembolism. Invest Radiol
18. Wang S, Unnikrishnan S, Herbst EB, et al. Ultrasound molecular imaging of inflammation in mouse abdominal aorta. Invest Radiol
19. Hyvelin JM, Tardy I, Bettinger T, et al. Ultrasound molecular imaging of transient acute myocardial ischemia with a clinically translatable P- and E-selectin targeted contrast agent: correlation with the expression of selectins. Invest Radiol
20. Tang MX, Mulvana H, Gauthier T, et al. Quantitative contrast-enhanced ultrasound imaging: a review of sources of variability. Interface Focus
21. Pitre-Champagnat S, Leguerney I, Bosq J, et al. Dynamic contrast-enhanced ultrasound parametric maps to evaluate intratumoral vascularization. Invest Radiol
22. Lassau N, Coiffier B, Faivre L, et al. Study of intrapatient variability and reproducibility of quantitative tumor perfusion parameters evaluated with dynamic contrast-enhanced ultrasonography. Invest Radiol
23. Desailly Y, Couture O, Fink M, et al. Sono-activated ultrasound localization microscopy. Appl Phys Lett
24. Christensen-Jeffries K, Browning RJ, Tang MX, et al. In vivo acoustic super-resolution and super-resolved velocity mapping using microbubbles. IEEE Trans Med Imaging
25. Hansen KB, Villagomez-Hoyos CA, Brasen JC, et al. Robust microbubble tracking for super resolution imaging in ultrasound
. IEEE International Ultrasonics Symposium Proceedings; 2016.
26. Desailly Y, Tissier AM, Correas JM, et al. Contrast enhanced ultrasound by real-time spatiotemporal filtering of ultrafast images. Phys Med Biol
27. Sboros V, Moran CM, Pye SD, et al. The behaviour of individual contrast agent microbubbles. Ultrasound Med Biol
28. Sboros V, Pye SD, MacDonald CA, et al. Absolute measurement of ultrasonic backscatter from single microbubbles. Ultrasound Med Biol
29. Sboros V, Pye SD, Anderson TA, et al. Acoustic Rayleigh scattering at individual micron-sized bubbles. Appl Phys Lett
30. Viessmann OM, Eckersley RJ, Christensen-Jeffries K, et al. Acoustic super-resolution with ultrasound and microbubbles. Phys Med Biol
31. Errico C, Pierre J, Pezet S, et al. Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging. Nature
32. Ackermann D, Schmitz G. Detection and tracking of multiple microbubbles in ultrasound B-Mode Images. IEEE Trans Ultrason Ferroelectr Freq Control
33. Opacic T, Dencks S, Theek B, et al. Motion model ultrasound localization microscopy for preclinical and clinical multiparametric tumor characterization. Nat Commun
34. Couture O, Hingot V, Heiles B, et al. Ultrasound localization microscopy and super-resolution: a state of the art. IEEE Trans Ultrason Ferroelectr Freq Control
35. Ackermann D, Schmitz G. Reconstruction of flow velocity inside vessels by tracking single microbubbles with an MCMC data association algorithm. Proc IEEE Int Ultrason Symp (IUS); 2013:627–630.
36. Foiret J, Zhang H, Ilovitsh T, et al. Ultrasound localization microscopy to image and assess microvasculature in a rat kidney. Sci Rep
37. Sboros V, Moran CM, Pye SD, et al. An in vitro study of a microbubble contrast agent using a clinical ultrasound imaging system. Phys Med Biol
38. Gabe IT, Gault JH, Ross J, et al. Measurement of instantaneous blood flow velocity and pressure in conscious man with a catheter-tip velocity probe. Circulation
39. Sugii Y, Okuda R, Okamoto K, et al. Velocity measurement of both red blood cells and plasma of in vitro blood flow using high-speed micro PIV technique. Measure Sci Tech
. 2005;16:1126. Available at: https://iopscience.iop.org/article/10.1088/0957-0233/16/5/011
. Accessed March 23, 2019.
40. Bar-Zion A, Solomon O, Tremblay-Darveau C, et al. Sparsity-based ultrasound super-resolution hemodynamic imaging. IEEE Trans Ultrason Ferroelectr Freq Control
41. van Sloun RJG, Solomon O, Bruce M, et al. Super-resolution ultrasound localization microscopy through deep learning. Ithaca, NY: Cornell University; 2018. Available at: https://arxiv.org/abs/1804.07661v2
. Accessed March 23, 2019.
42. Yang L, Qiu Z, Greenaway AH, et al. A new framework for particle detection in low-SNR fluorescence live-cell images and its application for improved particle tracking. IEEE Trans Biomed Eng
43. Wilson RS, Yang L, Dun A, et al. Automated single particle detection and tracking for large microscopy datasets. R Soc Open Sci
44. Yang L, Parton R, Ball G, et al. An adaptive non-local means filter for denoising live-cell images and improving particle detection. J Struct Biol
45. Butler MB, Thomas DH, Pye SD, et al. The acoustic response from individual attached and unattached rigid shelled microbubbles. Appl Phys Lett
46. Thomas DH, Looney P, Steel R, et al. Acoustic detection of microbubble resonance. Appl Phys Lett
47. Butler MB, Thomas DH, Silva N, et al. On the acoustic response of microbubbles in arteriole sized vessels. Appl Phys Lett
48. Thomas DH, Butler M, Anderson T, et al. The “quasi-stable” lipid shelled microbubble in response to consecutive ultrasound pulses. Appl Phys Lett
. 2012;101:071601. Available at: https://doi.org/10.1063/1.4746258
. Accessed March 26, 2019.
49. Thomas DH, Butler M, Pelekasis N, et al. The acoustic signature of decaying resonant phospholipid microbubbles. Phys Med Biol
50. Herbst EB, Unnikrishnan S, Wang S, et al. The use of acoustic radiation force decorrelation-weighted pulse inversion for enhanced ultrasound contrast imaging. Invest Radiol
51. Couture O, Besson B, Montaldo G, et al. Microbubble ultrasound super-localization imaging (MUSLI). 2011 IEEE International Ultrasonics Symposium, 18–21 October 2011
52. O'Reilly MA, Hynynen K. A super-resolution ultrasound method for brain vascular mapping. Med Phys
53. Sboros V, Averkiou M, Lampaskis M, et al. Imaging of the ovine corpus luteum microcirculation with contrast ultrasound. Ultrasound Med Biol
54. Pinker K, Moy L, Sutton EJ, et al. Diffusion-weighted imaging with apparent diffusion coefficient mapping for breast cancer detection as a stand-alone parameter: comparison with dynamic contrast-enhanced and multiparametric magnetic resonance imaging. Invest Radiol
55. Kapetas P, Clauser P, Woitek R, et al. Quantitative multiparametric breast ultrasound: application of contrast-enhanced ultrasound and elastography leads to an improved differentiation of benign and malignant lesions. Invest Radiol
56. Wang S, Herbst EB, Mauldin FW Jr., et al. Ultra-low-dose ultrasound molecular imaging for the detection of angiogenesis in a mouse murine tumor model: how little can we see? Invest Radiol
57. Lassau N, Bonastre J, Kind M, et al. Validation of dynamic contrast-enhanced ultrasound in predicting outcomes of antiangiogenic therapy for solid tumors: the French multicenter support for innovative and expensive techniques study. Invest Radiol
58. Hudson JM, Williams R, Karshafian R, et al. Quantifying vascular heterogeneity using microbubble disruption-replenishment kinetics in patients with renal cell cancer. Invest Radiol
59. Smeenge M, Tranquart F, Mannaerts CK, et al. First-in-human ultrasound molecular imaging with a VEGFR2-specific ultrasound molecular contrast agent (BR55) in prostate cancer
: a safety and feasibility pilot study. Invest Radiol
60. Huellner MW, Pauli C, Mattei A, et al. Assessment of prostate cancer
with dynamic contrast-enhanced computed tomography using an en bloc approach. Invest Radiol
61. Mischi M, Turco S, Lavini C, et al. Magnetic resonance dispersion imaging for localization of angiogenesis and cancer growth. Invest Radiol
62. Yang Y, Li X, Wang T, et al. Integrated optical coherence tomography, ultrasound and photoacoustic imaging for ovarian tissue characterization. Biomed Opt Express
63. Schindelin J, Rueden CT, Hiner MC, et al. The ImageJ ecosystem: an open platform for biomedical image analysis. Mol Reprod Dev
64. Schneider CA, Rasband WS, Eliceiri KW. NIH Image to ImageJ: 25 years of image analysis. Nat Methods
65. Grazul-Bilska AT, Borowicz PP, Reynolds LP, et al. Vascular perfusion with fluorescent labeled lectin to study ovarian functions. Acta Histochem
66. Boujelben A, Watson M, McDougall S, et al. Multimodality imaging and mathematical modelling of drug delivery to glioblastomas. Interface Focus
67. Martin K. Properties, limitations and artefacts of B-Mode images. In: Hoskins PR, Martin K, Thrush A, eds. Diagnostic Ultrasound: Physics and Equipment
. 2nd ed. Cambridge, United Kingdom: Cambridge University Press; 2010:64–74.
68. Ng A, Swanevelder J. Resolution in ultrasound imaging. Cont Educ Anaes Crit Care Pain
. 2011;11:186–192. Available at: https://doi.org/10.1093/bjaceaccp/mkr030
. Accessed March 26, 2019.
69. Diamantis K, Anderson T, Butler MB, et al. Resolving ultrasound contrast microbubbles using minimum variance beamforming. IEEE Trans Med Imaging
70. Kaplan D, Ma Q. On the statistical characteristics of log-compressed Rayleigh signals: theoretical formulation and experimental results. J Acoust Soc Am
71. Michailovich OV, Tannenbaum A. Despeckling of medical ultrasound images. IEEE Trans Ultrason Ferroelectr Freq Control
72. Perperidis A, Thomas D, Averkiou M, et al. Automatic dissociation between microvasculature and larger vessels for ultrasound contrast imaging. Conf Proc IEEE Eng Med Biol Soc
73. Song P, Trzasko JD, Manduca A, et al. Improved super-resolution ultrasound microvessel imaging with spatiotemporal nonlocal means filtering and bipartite graph-based microbubble tracking. IEEE Trans Ultrason Ferroelectr Freq Control
74. Yamanaka M, Smith NI, Fujita K. Introduction to super-resolution microscopy. Microscopy (Oxf)
75. Chenouard N, Smal I, de Chaumont F, et al. Objective comparison of particle tracking methods. Nat Methods
76. Montironi R, van der Kwast T, Boccon-Gibod L, et al. Handling and pathology reporting of radical prostatectomy specimens. Eur Urol
77. Bar-Zion A, Tremblay-Darveau C, Solomon O, et al. Fast vascular ultrasound imaging with enhanced spatial resolution and background rejection. IEEE Trans Med Imaging
78. Brawer MK. Quantitative microvessel density: a staging and prognostic marker for human prostatic carcinoma. Cancer
79. Russo G, Mischi M, Scheepens W, et al. Angiogenesis in prostate cancer
: onset, progression and imaging. BJU Int
80. Kerns MJJ, Darst MA, Olsen TG, et al. Shrinkage of cutaneous specimens: formalin or other factors involved? J Cutan Pathol
81. Leong AS, Leong TY. Newer developments in immunohistology. J Clin Pathol
82. Eisenthal A, Trejo L, Shtabsky A, et al. A novel assessment of the quality of immunohistostaining overcomes the limitations of current methods. Pathol Res Pract