Bubbly flows appear in various industrial applications and their measurement is an important topic. Specifically, the bubble size is an important parameter for determining the phase interfacial area and the transfer of mass, momentum and heat between phases. The bubble chord lengths can be measured by an intrusive technique such as resistance  or optical probe  . The bubble size can be obtained from bubble chord lengths using the bubble curvature  or statistical modelling  , but there is a bias towards larger bubbles, which have a higher probability of passing the probe than smaller bubbles. Bubble sizes can also be measured by the wire mesh tomography  with a limited resolution. The most common non-intrusive option is a high-speed camera (HSC). However, the HSC requires an optical access to each bubble, which is not possible in case of opaque pipes, liquids or large volumes of bubbly flow, when bubble overlapping occurs even for very low void fractions  . Techniques to detect bubble overlapping such as    are available, but they are still not adequate for measuring bubble sizes inside large volumes of flow. In this paper, we focus on ultrasound, which can be applied to these cases.
Ultrasound can measure bubble sizes using the resonance excitation  . However, this technique has a low accuracy, position resolution and is not suitable for measuring bubble size distributions. The Doppler technique  is suitable for velocity field measurement but does not detect the bubble size. The pulse echo technique measures the axial (along an ultrasonic beam) distance between the ultrasonic transducer (TDX) and a reflecting object. The radius of a single vapour bubble growing rapidly on a heated surface was measured by the pulse echo in  . It is adequate to measure only from one direction in cases with known position of the opposing bubble surface (e.g. plug bubbles in the horizontal pipe flow   ). Otherwise, measurements from two directions are necessary to get the bubble size  . Double array TDX measurement of bubble sizes at different vertical positions was proposed in  and realised in  . However, the pulse echo technique was limited to single bubbles in  -  . In the previous work  , the authors have shown that by adding a tracking technique the ultrasonic pulse echo can be applied to measure the interface positions of multiple bubbles simultaneously.
This paper aims to experimentally validate the performance of the ultrasonic pulse echo with tracking technique for the simultaneous measurement of multiple bubble sizes. The bubble size is represented by the horizontal length of the bubble, which was compared with the same parameter measured by a HSC. The paper shortly introduces the technique used in  with a single TDX and describes in detail how are data from two different TDXs combined to obtain bubble diameters. The experimental procedure is explained in detail including the video processing technique applied for analysing the HSC videos and the procedure used to compare HSC results with the pulse echo results.
2. Principle of the Measurement Technique
The technique applied the ultrasonic pulse echo with tracking technique  to measure the distance between a TDX and a bubble. The distances measured in each pulse repetition were tracked to form bubble surface trajectories. By combining the measurement from two opposing directions, two trajectories corresponding to opposing sides of the bubble were obtained. The length of the bubble along the axial direction of TDXs was calculated using the distance between TDXs x12 as
where x1 and x2 are the instantaneous TDX-bubble distances measured by TDX1 and TDX2 respectively (Figure 1). This section shortly introduces the pulse echo with tracking technique which detects bubble surface trajectories. Then, a calibration is explained. Finally, a technique combining trajectories from two TDXs to get bubble lengths is introduced. Figure 2 shows the flowchart of the whole technique.
2.1. Pulse Echo with Tracking Technique
The technique was fully described in  and is shortly introduced here. A TDX emitted ultrasonic pulses of N periods of basic frequency f0 in repetitions. The TDX detected waveforms with echo signals from bubbles and the waveform of each pulse was analysed to detect the transit time (delay) of bubble echoes in real time. The transit time T depends on the distance between the TDX and the bubble x and the speed of sound c as
and was estimated using a cross-correlation of the waveform with a reference signal (a measured reflection of the pulse from a free surface). This way, one bubble echo point was detected in the waveform of each pulse repetition for every bubble present at the sound field during that pulse’s emission.
Figure 1. The illustration of the technique. Each transducer (TDX) measures the distance between the bubble and itself and by subtracting these from the known distance between transducers x12, the bubble length d is obtained.
Figure 2. The flowchart of the measurement technique.
The detected bubble echo points represent positions of the bubble surface and were connected to trajectories using a tracking scheme described in detail in  . As a result, the technique provides a number of trajectories containing several points (t, x, M) where t stands for the measurement time (number of pulse in the sequence of pulses), and x with M are the position and echo signal strength respectively.
The relationship from Equation (2) is not exact because it does not include the time T0 when the sound wave travels through the TDX itself (e.g. the matching layer of the TDX). This results in a small error, which can usually be overlooked in ultrasonic applications. Because a high accuracy is necessary in this case, a calibration was conducted by measuring reflection from a steel prism placed in a known distance xc from the TDX surface. The measured echo delay Tc allows calculating the time T0 as
The value of T0 was subtracted from the transit time of each point. The speed of sound c was obtained from  using the measured temperature of the water.
2.3. Combining Trajectories from Opposing Transducers
The calculation of bubble length required that trajectories from opposing TDXs corresponding to the same bubble are paired together. The pairing used a distance function , which was low when there was a high probability that trajectories belong to the same bubble. In each step of the pairing algorithm, the trajectory pair i-j with the minimum distance function was paired and then both paired trajectories were removed from the list of considered trajectories.
The distance function was calculated as follows. First, positions measured by the TDX1 (distance from TDX1) were transformed to the distance from the TDX2 as
Then, for trajectory (with m points) measured by TDX1 and (with n points) measured by TDX2, the distance function was calculated as
using an assumed bubble length a, a scaling factor w and comparing the beginning, the end and the mean of each trajectory. The notation represents the k-th point of the i-th trajectory measured by the TDX1 and is the operation of averaging over all points of the i-th trajectory measured by the TDX1.
Figure 3 shows an example of trajectories for pairing. The correct pairs are 1-a, 2-b, 3-c, etc. To avoid incorrect pairings, such as the pair 2-a, which would result in a negative bubble volume, the was set to an unrealistically high value of 1000 s if to prioritise other possible connections of and overobviously incorrect pairs. An overlapping, such as the pair 6-c might occur if a shadow or interference affects echoes of two close bubbles (3-c and 6-f). After the bubble 6-f entered the measurement volume, the bubble 3-c was not detected by the TDX1 anymore (trajectory 3 ended) and the trajectory f was also affected. This can lead to a severe overestimation of bubble length. This problem is avoided by a proper choice of the assumed bubble length a and the factor w, which was set as 5 s/m because this reflected typical differences of 5 ms and 1 mm in time and position, respectively. The assumed bubble length was first set as 5 mm and then iteratively adjusted (while repeating the pairing of the whole dataset) to a median of the bubble length distribution. This approach becomes less effective for wide or non-uniform distribution of bubble lengths.
Despite the precautions described above, some pairs were incorrect. A common cause was a trajectory with no suitable trajectory from the opposite TDX nearby, a result of a bubble being detected only by one TDX or of the trajectories being just a noise. These lone trajectories were usually connected to other lone trajectories with a large difference in their measurement times and positions. These errors were identified by checking the length of a common measurement time interval for trajectory pairs. The pairs for which
were deleted. Finally, positions were interpolated on the common time interval to get new evenly sampled trajectories x’ with u points. The average horizontal bubble length d was calculated as the difference of average positions from both interpolated trajectories:
The pairs with d over 25 mm were considered errors and were deleted as well.
3. Experimental Validation of the Technique
The technique was validated using an experimental apparatus from Figure 4 and a concurrent HSC measurement as a reference. Air was supplied by an air compressor into an air feed pipe with a closed top end and valves regulated the air flow rate Fa. Air continued from the air feed pipe through a narrow nozzle of 0.4 mm in diameter into an annular section filled with water. Bubbles were formed at the nozzle outlet and moved upwards into an acrylic box filled with water and to the water surface. The annulus around the air feed pipe was 0.25 mm wide and water could be supplied to it at a valve regulated flow rate Fw. Because of the water flow in the annulus, an upwards drag force acted upon the bubbles forming at the nozzle outlet and speeded up the bubble detachment from the nozzle. As a result, the bubble volume was controlled using the water flow rate.
Figure 3. An example of trajectories for pairing.
Figure 4. Experimental apparatus. a) Front view as seen by the HSC1. b) Side view as seen by the HSC2.
A pulse generator/receiver JPR-10C-8CH3R (Japan Probe) with parameters set as in Table 1 generated pulses at a pulse repetition frequency fPRF = 1 kHz and with the second signal delayed to avoid interference of the ultrasonic pulses. Sound waves were generated and recorded by TDX1 (Imasonic 2M10) and TDX2 (Japan Probe B2Z10I). Waveforms filtered and amplified using a 200 kHz high-pass filter and 12 dB gain were digitised at the sampling speed of 200 MS/s by the APX-5040 (Aval Data) AD converter, which was deployed in a PC with Intel Core i5-4690 processor. The pulse echo signal was processed online in this PC, the tracking and pairing of trajectories were conducted later.
3.1. Reference Measurement by a High-Speed Camera
Two HSCs were used: Fastcam Mini AX50 (Photron) as HSC1 and Phantom Miro LC310 (Vision Research Inc.) as HSC2. The HSC1 was vertically aligned with the TDXs and it measured bubble lengths, which were compared with the ultrasound results for each detected bubble. The HSC2 was positioned just under the TDXs and it detected whether each bubble passes around or through the ultrasound beam and therefore whether it should be detected by the ultrasound or not.
Table 2 describes the HSC settings. Both cameras used backlight imaging. Bubbles cast shadows, which were clearly visible on the white background of the recorded videos. A grey-scale threshold was applied to obtain binary images from video frames in Matlab. Binary images (1 for bubble, 0 for no bubble) were filtered via morphological closing and bubble perimeters were detected using the Moore neighbour tracing algorithm with Jacob’s stopping criteria  .
Overlapping bubbles were detected and separated using techniques from  . The curvature was calculated from a filtered perimeter slope and concavities were detected as local curvature minima using a threshold of −0.1 px−1 (−0.6 mm−1) or −0.3 px−1 for large bubbles to avoid detecting small concavities, which might appear even on a single bubble. The bubble perimeter was separated at concavities and separated segments were fitted with ellipses  and the segments with overlapping ellipses were reconnected  . Finally, lines were drawn between reconnected segments to enclose separated bubbles. This can cause incorrect bubble shapes and volumes but improves the bubble length estimation as opposed to using the fitted ellipses such as in  .
Finally, bubbles were tracked using a PTV algorithm  and changes of the bubble centre position. Then, tracked bubbles from HSC1 were paired with tracked bubbles from HSC2. The distance function used for the pairing used differences in horizontal and vertical lengths, differences in vertical velocities and differences in the vertical position at the times when HSC1 trajectories begin (extrapolated from the HSC2 final position using the average of vertical velocities due to the difference in vertical positions of cameras resulting in about 15 ms gap between HSC2 and HSC1 trajectories).
Table 1. Parameters of generated pulses.
Table 2. Parameters of both applied HSCs.
3.2. Method for Comparing the Measurement Results
The HSC data for bubbles which did not pass the ultrasound field were deleted and the rest of the HSC data (pairs of trajectories from HSC1 and HSC2) were paired with the ultrasound data (pairs of trajectories from TDX1 and TDX2). The average horizontal bubble length measured by ultrasound was then compared with the same measured by the HSC1. The comparison was made for each bubble.
The decision of whether a bubble passes the ultrasound field or not was based on the prediction of its position perpendicular to the ultrasound beam (Figure 5). For each point of HSC1 trajectory, the position y was estimated using
where z1 was the vertical position measured by HSC1, y2l and z2l were horizontal and vertical positions of the last point l of the HSC2 trajectory. The was the standard deviation of all detected tilts which were calculated as
using the first and last horizontal and vertical coordinates of the HSC2 trajectory. The ultrasound beam boundary was evaluated numerically as in  (the shape roughly resembled a 10 mm wide cylinder) and the HSC1 trajectory points for which the complete interval (ymin, ymax) or the vertical position z1 were outside of this boundary were deleted. On average, 5% of HSC trajectories were deleted and trajectories were shortened to include only the part inside the ultrasonic beam.
Figure 5. An illustration of the decision whether a bubble passes the ultrasound (US) field. The figure is a superposition of 12 frames from HSC2 with a time step of 5.7 ms between them. The background is white, the TDX holder is black and bubbles are dark grey. Stationary noise (small dark circles) was emphasized by the superposition and was less obvious on individual images.
The pairing of HSC data with ultrasound data resembled the one in Section 2.3 with the distance function calculated as
where was the pair of ultrasonic trajectories with u points long interpolated trajectories and was a HSC trajectory pair with h points long HSC1 trajectory where positions and x2 were the maximum and minimum distance between bubble boundary and the TDX2 surface respectively. The scaling factor of w' = 5 s/m was used (same as in the Section 2.3). The pairs with were deleted as incorrectly connected.
4. The Performance Evaluation
The results of the experimental validation are presented together with the uncertainty analysis, which shows two types of error: common errors affecting all bubbles and uncommon errors which affect only some bubbles. An additional experiment using metal cylinders of known diameters was conducted and the measurement results were analysed to validate uncertainties for common errors. Then, uncommon errors were discussed including their dependency on the measurement conditions.
4.1. Comparison of Measured Bubble Lengths
The measurement was conducted for 6 measurement cases with different air and water flow rates and thus different bubble diameters, bubble number densities and vertical bubble velocities. The gas flow rate varied from 0.125 to 0.6 l/min, the water flow rate was less than 0.1 l/min, which was outside of the flow meter measurement range. The Morton number was 2.2 × 10−11, the Eotvos number was 0.3 - 5.6 and the Weber number 7.8 - 33.5. The void fraction was roughly estimated from the HSC data and ranged from 0.01% to 1.25% in different measurement cases. The measurement for each case was 29.105 s long and hundreds of bubble were detected and compared. An example of paired ultrasound trajectories for 3 simultaneously detected bubbles is shown in Figure 6 and a HSC1 image from the same time is shown in Figure 7. The data show that while one trajectory of TDX2 is shortened due to the corresponding bubble entering the shadow of another bubble, the presence of the bubble was still successfully detected, and its length was measured. The bubble length comparison for all detected bubbles is shown in Figure 8. The maximum number of bubbles detected simultaneously by the pulse echo was 6.
Figure 8 compares the measured bubble lengths. A histogram is used because there are over 1500 bubbles compared and the scatter plot does not clearly distinguish various high concentrations of points. Uncertainties of ±0.17 mm, 0.33 mm and 1.32 mm corresponded to a 50%, 68% and 95% percentile respectively and a slight overprediction was observed for small diameters. The difference distribution is not a normal distribution as can be understood from the difference of standard deviations (0.33 and 0.66 mm) implied by the 68% and 95% percentiles. The HSC1 resolution was 0.17 mm/px.
Figure 6. Measurement case 1: an example of 3 ultrasound trajectory pairs.
Figure 7. Measurement case 1: an example of a picture taken by the HSC1.
Figure 8. Comparison of horizontal bubble lengths for all bubbles. Scatter plot with void fractions estimated form the HSC data and 2D histogram.
4.2. Uncertainty Analysis
Various sources of error are listed in Table 3. Common errors affect all measured bubbles and are discussed further in Section 4.3. On the other hand, some uncommon errors affect only few measured bubbles. Uncommon errors and the ratio of bubbles affected by them are discussed in Section 4.4. This section explains all the listed sources of errors.
The uncertainty of the distance x12 between TDXs includes the measurement uncertainty of x12 and changes due to thermal expansion, vibration, etc. The uncertainty in the speed of sound calculation includes the temperature measurement uncertainty (±1 K) and variations in both time and space (assumed as ±0.5 K in this experiment). These, together with the measurement uncertainty of T0 were systematic errors (considering a short measurement period during which neither x12 nor c changes significantly).
The uncertainty of the transit time T calculation as the cross-correlation maximising time lag depends on the sampling speed and on the noisiness of the data. The cross-correlation function would ideally become the autocorrelation function of the reference signal (Figure 9), which resembles the cosine function. Assuming the width at the half maximum 1/(6f0) as the uncertainty of T, the position estimation uncertainty becomes 0.06 mm for one TDX and the bubble length uncertainty becomes 0.09 mm. A statistical test could provide a more accurate uncertainty estimation.
Table 3. Sources of measurement error in the pulse echo measurement of bubble length.
Figure 9. Autocorrelation function of the reference signal.
Figure 9 illustrates another error of the transit time T measurement as well. In a noisy signal, a secondary peak at 1/f0 or a tertiary peak at 2/f0 can be mistaken for the main peak. This issue was the main topic of  . In this work, this error was treated by a mode filter described in  . However, the mode filter fails for too noisy signals and errors of 0.36 mm or its multiples appear. A strong noise can occur due to an interference with multiple reflections of the pulse previously emitted by the opposing TDX or due to a shadow cast by another bubble. Finally, an incorrect pairing of trajectories can cause severe over prediction of the bubble length similar to the overlapping error of HSC video processing. Unlike previously discussed errors, these uncommon errors appeared only in few bubbles and most bubbles were not affected by these errors.
4.3. Common Errors
The errors which affect all detected bubbles add up to the total measurement uncertainty of 0.31 mm. This estimation was validated by measuring the diameter of several metal cylinders, which were manually held at the centre of the measurement volume. The diameters were measured by callipers (accuracy of ±0.05 mm) as well and Table 4 compares the results. Differences were smaller than 0.1 mm, which agrees with 0.1 mm reported previously in  for falling metal balls and a similar technique.
The bubble length results were analysed as well. Considering the HSC1 resolution of 0.17 mm/px, the expected measurement uncertainty of bubble length differences is 0.35 mm. The uncertainty of 0.33 mm corresponding to the 68% percentile is smaller and thus it can be concluded that the pulse echo measurement technique was successfully validated by the HSC measurement.
4.4. Uncommon Errors
Uncommon errors affect only some bubbles but are larger than common errors. Because the uncommon errors are at least 0.36 mm (Section 4.2), they can be distinguished from smaller common errors. Assuming only common errors, the percentage of bubbles with difference larger than 2 × 0.33 mm should be 5%, but the measured value was 15%. This indicates that 10% of all detected bubble lengths were affected by large uncommon errors.
Because uncommon errors are often the consequence of measuring multiple bubbles simultaneously, the effect of void fraction on these errors was analysed (Figure 10). The average difference between pulse echo and HSC measured lengths was smaller than 0.2 mm for each measurement case, which can be explained by common systematic errors (Table 3). The uncertainties corresponding to the 68% and 95% percentile grow with the void fraction and the growth is more obvious for the 95% data. This is explained by the ratio of large uncommon errors (calculated the same way as above), which was higher for higher void fractions. While measuring only single bubbles (such as in  -  ) or single metal cylinders (Section 4.3), it is easier to avoid interference and reduce the noisiness of the signal, and thus to avoid the uncommon errors. With increasing the void fraction, problems with interference, shadowing etc. will lead to increasing rate of uncommon error occurrence in the measurement.
Finally, an over prediction of few small bubble lengths by the pulse echo (Figure 8) is discussed. For each case, the average difference was subtracted from all the differences to get rid of common systematic errors and the standard deviation σ was calculated for each measurement case. 16% of all data should be larger than σ but the ratio increases for bubbles smaller than 3 mm (Figure 11). This shows a systematic error which affects 6.5% of bubbles in the measurement case with the smallest bubbles and fewer for larger bubbles. This error was caused by an incorrect pairing of trajectories discussed in Section 2.3 or by lower signal to noise ratio due to smaller amplitude of reflection from small bubbles  .
Figure 10. Analysis of void fraction effect on the measured differences for all measurement cases.
Figure 11. The ratio of bubble length differences larger than the standard deviation of the data for each measurement case. The ratio should theoretically be 16% (dotted line).
Table 4. Comparison of metal cylinder diameters measured by pulse echo and by callipers.
The paper shows that the ultrasonic pulse echo with tracking technique can measure bubble diameters in the range of 1 to 10 mm. Two types of errors were discussed. Common errors affect all bubbles with an uncertainty of 0.31 mm as confirmed by an additional experiment and a comparison with a reference HSC measurement. Uncommon errors cause a higher uncertainty but affect only a small ratio of bubbles. The comparison with HSC data showed that the ratio of affected bubbles increased with a void fraction from about 3% to about 18% in the currently tested range of bubbly flow conditions.
The main limitation of the technique is that it measures the bubble length only in one direction, which does not represent bubble sizes perfectly. An additional sensor pair would allow measuring the bubble length in another direction, but with worse time resolution and much smaller measurement volume (only the cross-section of two ultrasound fields). An equivalent diameter of the bubble, which represents bubbles more adequately, can be estimated by assuming an oblate shape of bubbles and using a correlation for the aspect ratio (e.g. from  ). However, such correlations have a high uncertainty. In addition, the procedure requires an assumption that the horizontal bubble length is the major axis of the bubble, which can be inaccurate in case of strong wobbling. All these issues obviously become less relevant in case of more spherical bubbles, such as for smaller bubbles (lower Eotvos number) or gas-liquid combination with a higher Morton number. In these cases, the one-dimensional bubble length represents bubbles adequately.
The greatest advantages of this technique are the ability to conduct a non-intrusive measurement of bubble sizes even in opaque environments, and the ability to immerse the sensors and thus measure bubbly flows in large geometries (not only in narrow channels) and the relatively low cost. These advantages make the ultrasonic pulse echo with tracking technique an attractive option for measuring bubble sizes in a two-phase flow.
 Jones, O.C. and Delhaye, J.M. (1976) Transient and Statistical Measurement Techniques for Two-Phase Flows: A Critical Review. International Journal of Multiphase Flow, 3, 89-116.
 Guet, S., Fortunati, R.V., Mudde, R.F. and Ooms, G. (2003) Bubble Velocity and Size Measurement with a Four-Point Optical Fiber Probe. Particle & Particle Systems Characterization, 20, 219-230.
 Clark, N.N. and Turton, R. (1988) Chord Length Distributions Related to Bubble Size Distributions in Multiphase Flows. International Journal of Multiphase Flow, 14, 413-424.
 Prasser, H.M., Scholz, D. and Zippe, C. (2001) Bubble Size Measurement Using Wire-Mesh Sensor. Journal of Flow Measurement and Instrumentation, 12, 299-312.
 Lacuona, A., Sosa, P.A., Rodriguez, P.A. and Zequeira, R.I. (2000) Volumetric Characterisation of Dispersed Two-Phase Flows by Digital Image Analysis. Measurement Science and Technology, 11, 1152-1161.
 Shen, L., Song, X., Iguchi, M. and Yamamoto, F. (2000) A Method for Recognizing Particles in Overlapped Particle Images. Pattern Recognition Letters, 21, 21-30.
 Honkanen, M., Saarenrinne, P., Stoor, T. and Niinimaki, J. (2005) Recognition of Highly Overlapping Ellipse-Like Bubble Images. Measurement Science and Technology, 16, 1760-1770.
 Fu, Y. and Liu, Y. (2016) Development of a Robust Image Processing Technique for Bubbly Flow Measurement in a Narrow Rectangular Channel. International Journal of Multiphase Flow, 84, 217-228.
 Fairbank, W.M. and Scully, M.O. (1977) A New Noninvasive Technique for Cardiac Pressure Measurement: Resonant Scattering of Ultrasound from Bubbles. IEEE Transactions on Biomedical Engineering, BME-24, 107-110.
 Murakawa, H., Kikura, H. and Aritomi, M. (2008) Application of Ultrasonic Multi-Wave Method for Two-Phase Bubbly and Slug Flows. Flow Measurement and Instrumentation, 19, 205-213.
 Duffey, R.B. and Hall, R.S. (1969) An Ultrasonic Technique for Measuring the Transient Movement of a Liquid-Vapour Interface. Journal of Physics E: Scientific Instruments, 2, 193-194.
 Matikainen, L., Irons, G.A., Morala, E. and Chang, J.S. (1986) Ultrasonic System for the Detection of Transient Liquid/Gas Interfaces Using the Pulse-Echo Technique. Review of Scientific Instruments, 57, 1661-1666.
 Chang, J.S. and Morala, E.C. (1990) Determination of Two-Phase Interfacial Areas by an Ultrasonic Technique. Nuclear Engineering and Design, 122, 143-156.
 Poole, D.H. (1970) An Ultrasonic Technique for Measuring Transient Movements of a Phase Boundary in a Liquid. Journal of Physics E: Scientific Instruments, 3, 726-728.
 Andruszkiewicz, A., Eckert, K., Eckert, S. and Odenbach, S. (2013) Gas Bubble Detection in Liquid Metals by Means of the Ultrasound Transit-Time-Technique. The European Physical Journal Special Topics, 220, 53-62.
 Richter, T., Eckert, K., Yang, X. and Odenbach, S. (2015) Measuring the Diameter of Rising Gas Bubbles by Means of the Ultrasound Transit Time Technique. Nuclear Engineering and Design, 291, 64-70.
 Povolny, A., Kikura, H. and Ihara, T. (2018) Ultrasound Pulse-Echo Coupled with a Tracking Technique for Simultaneous Measurement of Multiple Bubbles. Sensors, 18, 1327.
 Reddy, P.R., Amarnadh, V. and Bhaskar, M. (2012) Evaluation of Stopping Criterion in Contour Tracing Algorithms. International Journal of Computer Science and Information Technology, 3, 3888-3894.
 Fitzgibbon, A., Pilu, M. and Fisher, R.B. (1999) Direct Least Square Fitting of Ellipses. IEEE Transactions on Pattern Analysis and Machine Intelligence, 21, 476-480.
 Chang, T.P.K., Watson, A.T. and Tatterson, G.B. (1985) Image Processing of Tracer Particle Motions as Applied to Mixing and Turbulent Flow I. Chemical Engineering Science, 40, 269-275.
 Wada, S., Furuichi, N. and Shimada, T. (2017) Application of Partial Inversion Pulse to Ultrasonic Time-Domain Correlation Method to Measure the Flow Rate in a Pipe. Measurement Science and Technology, 28, Article ID: 115302.
 Besagni, G. and Inzoli, F. (2016) Bubble Size Distributions and Shapes in Annular Gap Bubble Column. Experimental Thermal and Fluid Science, 74, 27-48.