Bubble Tracking Algorithm

Below is an excerpt from my master's thesis Lateral Fluid Motion in Nucleate Boiling that explains how the velocity field of the bubbles was resolved by processing the high-speed camera frames to measure that change in position of bubbles over the known time between frames.

The full thesis can be downloaded from Oregon State University's website for free here: download full thesis. This section begins on the PDF document page 97 (thesis page 77)

6.1. Determination of Bubble Kinematics

In order to obtain quantitative results for parameters like velocity, acceleration, and bubble departure diameter, an image-processing algorithm was developed in Matlab®. The primary focus of image processing described in this section was to resolve the individual bubble velocities through bubble tracking. Briefly, bubble-tracking velocimetry entails finding the displacement of each bubble from one frame (A) to the next (B) and then dividing by the elapsed time to resolve the velocity. The bubble was first located within an image and positively identified. Once positively identified, each bubble in the first frame was paired with the corresponding bubble in the second frame. Any unmatched pairs were deleted and resulting bubble pairs were converted to individual vectors and displayed on a full field of view image of the resolved bubbles. The following paragraphs describe the image processing steps in more detail.

After loading the image, the image histogram was stretched linearly over the entire grayscale spectrum such that the maximum amount of contrast was achieved without washing out any pixels (Fig. 6.1). Additionally, this ensured the images used in analysis always had a similar histogram, improving repeatability.

The first step in resolving the bubbles was to apply edge detection filters. Edge detection filters are used to restrict the algorithm to tracking only bubbles within a two-dimensional plane (i.e. the plane of focus). Since the K2/SC lens has a very narrow depth of field (as low as 0.02 mm at full aperture), the bubbles located outside the depth of field appear greatly blurred, and are not resolved by the edge detection filters (depending on the sensitivity set by the user). Three different filters – gradient, Laplace, and Sobel – were applied to the image and the results of all filters were merged in the end. Using multiple filters allowed certain bubbles to be resolved that individual filters may have missed. The perimeters of bubbles still attached to the surface could not be resolved as there was no distinct edge between the bubble and the surface to resolve; only departed bubbles could therefore be resolved.

Figure 6.2 illustrates the steps of processing an image using the Sobel filter. Similar processing steps were performed for the two other filters and the results were combined prior to tagging the bubbles for tracking. Upon applying the edge detection filter (Fig. 6.2b), the next step was to binarize the image (Fig. 6.2c). All pixels with a value above the user-defined binarization threshold were set to a value “1” and the remaining pixels were set to a value of “0”. The optimal threshold value was determined iteratively for each of the three filters and kept constant for all experiments to allow for meaningful comparisons. Henceforth, the term “region” will be used to refer to any continuous group of binary pixels with a value “1”. Following binarization, the image was further processed to resolve as many regions as possible to represent the location and shape of bubbles in the image. Since the edge detection filters only resolved bubble perimeters, the shape of the resulting region is in the binarized image was usually a ring. For this reason all holes in the regions are filled to obtain regions representative of a complete bubble (Fig. 6.2d).

It is possible that some bubble edges were not fully resolved (i.e. there is a gap in the edge-detected ring). To retain these regions, a convex hull process was applied as shown in Fig. 6.2e. The convex hull refers to the smallest envelope containing the region with continuous convex geometry, and can also be thought of as the outline that a rubber band would trace if it were pulled around the region. If the set of points representing the convex hull bridges a gap that exceeds a threshold (set as a fraction of the region’s major diameter), the original bubble edge is retained. However, if the threshold is not exceeded, the region is replaced with its convex hull, effectively enclosing the incomplete edge of the bubble. Figure 6.2e shows that the smaller bubbles are resolved using the convex hull whereas the larger bubble to the lower left hand corner remains unresolved. The next image-processing step was deletion of regions that were artifacts resulting from the background or from noise. These regions were filtered simply based on an area threshold, specified by the user (Fig. 6.2f). Any regions connected to the image border were also removed at this stage, as there would be no way of resolving partial bubbles that move in or out of the frame in subsequent images.

After the regions were resolved, they were processed to retain only those regions that represented the location and shape of bubbles in the image. Two criteria had to be met for a region to “pass” as a bubble: (1) region solidity (area/convex area) > threshold – This criterion was applied to eliminate regions resolved without continuous convex edges, and (2) region eccentricity < threshold – This criterion was applied to eliminate regions with high eccentricity. Figure 6.3 shows which regions were eliminated for not meeting the solidity and eccentricity criteria, labeled red and green respectively. Such regions are unlikely to represent the outline of a bubble, which is usually spherical in shape. All threshold values were set by the user and could be adjusted to resolve bubbles whose shape would otherwise not permit classification as a bubble.

Following identification of bubbles in each image, bubble tracking was performed. The first step in bubble tracking was determining the locations of each region’s centroid. For any region with equivalent diameter equal to or greater than 0.1 mm, the centroid was computed as the centroid of the pixels in the region. For small regions with equivalent diameters less than 0.1 mm) the centroid location was weighted by the intensity of the corresponding pixels in the original image. Intensity weighting the pixels within the region ensured that any misrepresentation of the dot due to poor focus or incorrect binary threshold settings would have minimal effect on the accuracy of the centroid computation. Due to the fairly low quantity of bubbles in each image, it was not necessary to use advanced algorithms for determining bubble pairs between images. A nearest-neighbor approach was used to match bubbles in image A to those in image B. Figure 6.4 shows the bubbles pairs identified in image A and image B. For each region in image A, the distance between its centroid and each region’s centroid in image B was calculated. The pair with the shortest centroid-to-centroid distance was logged as a matched pair; their respective indexes in the region-numbering scheme were logged in a separate table. Each matched pair was filtered for violation of the area ratio threshold (area(A)/area(B) < threshold). This filter ensured that only bubbles with similar area would be logged as pairs. Those regions that failed this test were tagged orange and subsequently removed. If bubbles were tracked near the cavity in subcooled conditions, the area ratio threshold was set high so that the collapse of a bubble would not exclude it from analysis. If only the bubbles in the plume were of interest, then this ratio was set low, as bubble collapse was not observed to occur in the plume.

It is possible that two regions in image A are matched to a single region in image B. This often occurred because the second region in image B was not resolved, likely because the bubble moved out of the depth of field. In this case, the pair with the shortest centroid-to-centroid distance was retained, as this pair most likely represented the slight shift in position the bubble has made between image A and image B. Figure 6.5 shows a typical vector plot that results from identification of bubbles in images A and B. The length of the vectors is proportional to the square root of the bubble speed for better viewing of all vectors in the speed range.

Because the shape of the individual regions was rarely circular, the diameter of each bubble is computed an equivalent diameter for a circle with area equal to that of the region. 

Equation 6.2 provides a good estimate of the diameter for the purpose of comparing the data to existing models of bubble departure diameter. The bubble volume was computed assuming the region represents a cross-section of an oblate spheroid (i.e. the minor axis is the axis of symmetry for the volume):

The length of the major and minor axis used in Eq. 6.2 was computed automatically by Matlab® and given in terms of pixels. This was deemed the most accurate method to estimate bubble volume knowing only the shape of the cross section, and was considered quite accurate for bubbles appearing circular or elliptical in shape. However, for greatly deformed bubbles this is likely a meaningless quantity, as there is no axis of symmetry in such a case.

Additional error is expected from the bubble diameter calculations, as there was no means of determining the accuracy of the equivalent diameter estimates. Figure 6.6 shows a cropped portion of an image with five regions as resolved by the gradient, Laplace, and Sobel filters (red, green, and blue, respectively).

The gradient and Sobel filters often agreed quite well, whereas the Laplace filter usually resolved smaller regions. The parameters that controlled the size of the resolved region were the binarization threshold and the filter operator size; lower binary thresholds allowed fainter pixels to be included in the region and smaller filter operator dimensions increased the sensitivity of the filter, meaning only sharp gradients would pass as edges. The most effective binarization threshold and filter operator size were determined iteratively for each filter, and kept constant for all experiments. However, the ambiguity as to which filter is more accurate indicates that there is likely some error associated with estimating bubble size using these filters. Since the final region size was a combined result of all three filters, in most cases this simply meant the largest of the three. However, as there was no means of comparing the resolved data to the true values, the uncertainty for this error could not be computed. A flowchart describing the algorithm along with the purpose of the individual functions used is provided in Appendix C.