3D Reconstruction of a Single Bubble in Transparent Media Using Three Orthographic Digital Images

: This work proposed a method to reconstruct the 3D bubble shape in a transparent medium utilizing the three orthographic digital images. The bubble was divided into several ellipse slices. The azimuth angle and projection parameters were extracted from the top view image, while the formulas for dimensionless semi-axes were derived according to the geometric projection relationship. The elliptical axes of each layer were calculated by substituting the projection width into the formulas. All layers of slices were stacked to form the 3D bubble shape. Reconstruction accuracy was evaluated with spheres, ellipsoids, and inverted teardrops. The results show that the position contributes greatly to the reconstruction accuracy of the bubbles with serious horizontal deformation. The method in Bian et al. (2013) is sensitive to both horizontal and vertical deformations. The vertical deformation has little inﬂuence on the method in Fujiwara et al. (2004), whereas the horizontal deformation greatly impacts its accuracy. The method in this paper is negligibly a ﬀ ected by vertical deformation, but it does better in reconstructing single bubbles with large horizontal deformation. The azimuth angle a ﬀ ects the accuracy of the methods in Bian et al. (2013) and Fujiwara et al. (2004) more than the method in this paper.


Introduction
Studying the hydrodynamics of a single bubble growth and detachment in motionless fluids always attracts researchers [1,2]. The literature of a single bubble generated in Newtonian fluids has been studied for decades and this well-known classic work was published, [3]. Despite the importance of this problem, many applications are dealing with bubbles in non-Newtonian fluids, and bubbles in such fluids have not been thoroughly studied. Recently, some studies focused on the dependence of the single bubble shape on their dynamics such as velocity discontinuity, as well as on the contact angle, viscosity, and surface tension [4,5]. To this end, the 3D bubble shape reconstruction is vital for the further analysis of hydrodynamic interactions among a single bubble and surrounding liquid.
There exists a significant body of scientific research discussing the calculation of bubble size and its formation. Early researchers used the balance formula of buoyancy and surface tension to calculate the volume of a single bubble [6]. Then, the Rayleigh-Plesset equation [7,8] and Young-Laplace equation [9,10] were utilized to predict the formation of bubbles. Many numerical simulation techniques, including the boundary integral method [11], the finite element analysis method [12] have also been applied to forecasting bubble formation. In general, the results of theoretical calculations and numerical simulations need to be verified by accurate experimental data. Therefore, a method which could precisely extract 3D bubble formation from orthographic images can be employed to validate the rationality of various theoretical calculation methods and numerical simulation techniques. VCXG-13M), all connected to a computer, were utilized simultaneously for the front view, side view, and top view. Each camera had a matrix of 1024 × 1280 pixels and the size of the pixel was 0.0048 mm. The diameter of sparger was 1 mm. The distance between the front view camera and the center of the box was 463.0 mm, while the distance between the side view camera and the center of the box was 439.0 mm. Since the distance were much larger than experimental bubble diameter (less than 20 mm). Therefore, the distortion caused by light refraction between transparent soil and acrylic as well as the light refraction between acrylic and air was negligible. Two LED panels were placed behind the acrylic box to ensure sufficient illumination, which could reduce the significant errors in the bubble size measurement. Since the wall influence, which caused bubble asymmetry was uniform in cuboid acrylic box, the front view was similar to the side view and the top view resembled being circular.

Experimental Set-Up
The schematic of the image acquisition device is illustrated in Figure 1. Laponite RD transparent soil [28] was used as the supporting medium to maintain the suspension of a single bubble in the center of an acrylic box, the size of which was 100 × 100 × 150 mm 3 . Three digital cameras (Baumer VCXG-13M), all connected to a computer, were utilized simultaneously for the front view, side view, and top view. Each camera had a matrix of 1024 × 1280 pixels and the size of the pixel was 0.0048 mm. The diameter of sparger was 1 mm. The distance between the front view camera and the center of the box was 463.0 mm, while the distance between the side view camera and the center of the box was 439.0 mm. Since the distance were much larger than experimental bubble diameter (less than 20 mm). Therefore, the distortion caused by light refraction between transparent soil and acrylic as well as the light refraction between acrylic and air was negligible. Two LED panels were placed behind the acrylic box to ensure sufficient illumination, which could reduce the significant errors in the bubble size measurement. Since the wall influence, which caused bubble asymmetry was uniform in cuboid acrylic box, the front view was similar to the side view and the top view resembled being circular.

Image Processing and Camera Calibration
Images were processed by using MATLAB R2015a to obtain the bubble projection contours during the growth. The key result is shown in Figure 2. Figure 2a is one of the images taken by left view bubble in our experiment in Section 4.2. Figure 2b shows the image after dividing by a background image, which was without bubble, to enhance the contract and then binarizating with a threshold on the gray levels. From Figure 2c, we can notice that the holes are filled and noises are eliminated effectively. Bubble projection contours were detected by edge extraction method, as shown in Figure 2d.

Image Processing and Camera Calibration
Images were processed by using MATLAB R2015a to obtain the bubble projection contours during the growth. The key result is shown in Figure 2. Figure 2a is one of the images taken by left view bubble in our experiment in Section 4.2. Figure 2b shows the image after dividing by a background image, which was without bubble, to enhance the contract and then binarizating with a threshold on the gray levels. From Figure 2c, we can notice that the holes are filled and noises are eliminated effectively. Bubble projection contours were detected by edge extraction method, as shown in Figure 2d. After image processing, a calibrated chessboard was placed inside the acrylic box to determine the intrinsic parameters of cameras. According to the perspective of Zhang et al. and the application of their methods [29,30], the front and side camera parameters could be calculated, whose results were shown in Tables 1 and 2, separately. According to the ideal pinhole model, the magnification factors are determined then for front camera (0.177 mm pixel −1 ) and for side camera (0.171 mm pixel −1 ). Because the top view camera was only used to extract the azimuth angle, which cannot be effected by magnifying, therefore it did not need to be calibrated. After image processing, a calibrated chessboard was placed inside the acrylic box to determine the intrinsic parameters of cameras. According to the perspective of Zhang et al. and the application Appl. Sci. 2020, 10, 5803 4 of 18 of their methods [29,30], the front and side camera parameters could be calculated, whose results were shown in Tables 1 and 2, separately. According to the ideal pinhole model, the magnification factors are determined then for front camera (0.177 mm pixel −1 ) and for side camera (0.171 mm pixel −1 ). Because the top view camera was only used to extract the azimuth angle, which cannot be effected by magnifying, therefore it did not need to be calibrated. where, f ij is the calibrated focal length in pixels, subscript character i = 1,2 representing horizontal direction and vertical direction of pictures respectively while subscript character j = 1,2 representing front camera and side camera, respectively. γ is the skew coefficient defining the angle between the x and y pixel axes, k c is the image distortion coefficients (radial and tangential distortions).
The calibrated focal length was slightly different from the data provided from manufacturer that the camera was equipped with lens of focal distance 12 mm. There are probably two main reasons causing this unexpectedly disparity. On the one hand, ideal pinhole image model was adopted in camera calibration process, assuming focal length was definite. However, the real cameras had much more complex construction. On the other hand, different locations and positions of the calibration chessboard can result in variant calibration results. Thus, calibrating results change with distances and postures. Therefore, in order to receive a more satisfying calibration result, each step of calibration should be done exactly and carefully. Besides, the high quality of calibration chessboard, appropriate illumination conditions, and suitable camera postures should also be taken into consideration.

Measurement Accuracy
All steps of visualization produced uncertainties. First, high precision monox calibration chessboard was used to improve the measurement accuracy, the accuracy reached ±0.01 mm. Then, in order to eliminate the error produced by limited calibration pictures, this paper considers 30 different calibration postures as plentiful enough. We can easily know that one pixel created an measurement error of 0.0048 mm. Thus, the errors of hole filling as well as edge extraction could be measured and each of them exerted a negligible 0.1% on the measurement bubble length, as background division and binarization. Further, both front and side cameras induced uncertainty of local length through camera calibration process, which was [37.48217, 40.70242] and [56.57221, 64.13330] in pixels. Hence, the calibration errors in x and y pixel axes cold be calculated, respectively. In this way, the uncertainty on the measurement of bubble diameter was accumulated. According to the propagation law of errors: where, m d is the objective mean square error of bubble diameter, D j is the distance between the cameras and the center of the box, and was measured by standard ruler. Therefore, the mean square error m D was ±1 mm on the basis of ruler precision, m f is the calibration mean square errors, px is the length of a pixel, m c is the precision of calibration chessboard. Thus, after all the above operations, the measurement uncertainty on bubble diameter during camera calibration and image processing was estimated to be ±0.46 mm. The mean square error of spherical bubble volume can be deduced roughly as following: where, m V is the objective mean square error of bubble volume, of which the measurement is the uncertainty evolution of bubble geometrical characteristics. Then the relative mean square error (K) could be written as: where, V b is the volume of spherical bubbles. It can be seen in Figure 3 that K was the decreasing function of bubble diameter. With the bubble growth, the K met a relative stable level and the measurement accuracy increased. Pleased results (given the expected threshold K3%) could be met in the interval [4, +∞], which was quite acceptable in our experiments, where most of the bubble diameters were in this interval.
where, md is the objective mean square error of bubble diameter, Dj is the distance between the cameras and the center of the box, and was measured by standard ruler. Therefore, the mean square error mD was 1  mm on the basis of ruler precision, mf is the calibration mean square errors, px is the length of a pixel, mc is the precision of calibration chessboard. Thus, after all the above operations, the measurement uncertainty on bubble diameter during camera calibration and image processing was estimated to be 46 . 0  mm. The mean square error of spherical bubble volume can be deduced roughly as following: where, mV is the objective mean square error of bubble volume, of which the measurement is the uncertainty evolution of bubble geometrical characteristics. Then the relative mean square error (K) could be written as: where, Vb is the volume of spherical bubbles. It can be seen in Figure 3 that K was the decreasing function of bubble diameter. With the bubble growth, the K met a relative stable level and the measurement accuracy increased. Pleased results (given the expected threshold K≦3%) could be met in the interval [4, +∞], which was quite acceptable in our experiments, where most of the bubble diameters were in this interval.

Modified Stacking Ellipse Method
The method in this paper proposed here included four steps as illustrated in Figure 4. The first step was to extract the projected contour of the bubble in three orthographic images, and the second step corrected the posture of the projected contour to make the c-axis (or z-axis) in the front view and the side view vertical or horizontal; the third step was also supposed to calculate the dimensionless major semi-axis (a0) and minor semi-axis (b0) of the horizontal elliptical cross-section by utilizing the projection relationship and the geometric information; the fourth step divided the bubble into n layers along the height direction (c-axis) and calculated the major semi-axis (ai) and the minor

Modified Stacking Ellipse Method
The method in this paper proposed here included four steps as illustrated in Figure 4. The first step was to extract the projected contour of the bubble in three orthographic images, and the second step corrected the posture of the projected contour to make the c-axis (or z-axis) in the front view and the side view vertical or horizontal; the third step was also supposed to calculate the dimensionless major semi-axis (a 0 ) and minor semi-axis (b 0 ) of the horizontal elliptical cross-section by utilizing the projection relationship and the geometric information; the fourth step divided the bubble into n layers along the height direction (c-axis) and calculated the major semi-axis (a i ) and the minor semi-axis (b i ) of each horizontal ellipse slice. Finally, all the ellipse slices were stacked to form a 3D bubble shape.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 18 semi-axis (bi) of each horizontal ellipse slice. Finally, all the ellipse slices were stacked to form a 3D bubble shape.

Bubble Posture Correction
Bubble posture correction was the pre-work of the method proposed in this paper. There were two major purposes to carry out this step. On the one hand, based on the rotation conditions and geometric information of front contour and side contour around x-axis and y-axis, the posture correction improved the measurement accuracy of the azimuth angle around z-axis. On the other hand, the posture correction, which used the spatial rotation angles of bubble to correct their projection contours, properly reduced the measurement error causing by rotation. Figure 5 illustrates the schematic diagram of posture correction. Since the image processing steps of the front view are similar to those of the side view, as an example, the following statements are used to explain the steps of posture correction for the front view. In the following expressions, F, S, and T indicate the front view, side view, and top view, respectively.
The implementation details of bubble posture correction were shown in Figure 5. We considered images of sample No. 7 (As shown in the bubble examples in Section 4.1) to be the example. Firstly, the projected contours were extracted using image processing software MATLAB R2015a. Axis A was assumed to be the shorter axis, which is perpendicular to the longer axis, i.e., axis B. Moreover, the angle between axis A and the vertical direction, i.e., θF, was considered to be in the range of 0-90° (0° ≤ θF ≤ 90°). When θF is smaller than 45°, the projected contour of the top view of the bubble is mainly affected by the length of axis B, so axis A is rotated to the vertical direction; on the other hand, when θF is greater than 45°, the projected contour of the top view of the bubble is mostly affected by axis A, so axis A is rotated to the horizontal direction. The rotation formula is given by: where, (xF, zF) and ( F x , F z ) represent the coordinates of the bubble projection contour before rotation and after rotation, respectively; θF is the angle of rotation. Then, the magnification factor of the front view is calculated by:

Bubble Posture Correction
Bubble posture correction was the pre-work of the method proposed in this paper. There were two major purposes to carry out this step. On the one hand, based on the rotation conditions and geometric information of front contour and side contour around x-axis and y-axis, the posture correction improved the measurement accuracy of the azimuth angle around z-axis. On the other hand, the posture correction, which used the spatial rotation angles of bubble to correct their projection contours, properly reduced the measurement error causing by rotation. Figure 5 illustrates the schematic diagram of posture correction. Since the image processing steps of the front view are similar to those of the side view, as an example, the following statements are used to explain the steps of posture correction for the front view. In the following expressions, F, S, and T indicate the front view, side view, and top view, respectively.
The implementation details of bubble posture correction were shown in Figure 5. We considered images of sample No. 7 (As shown in the bubble examples in Section 4.1) to be the example. Firstly, the projected contours were extracted using image processing software MATLAB R2015a. Axis A was assumed to be the shorter axis, which is perpendicular to the longer axis, i.e., axis B. Moreover, the angle between axis A and the vertical direction, i.e., θ F , was considered to be in the range of 0-90 • (0 • ≤ θ F ≤ 90 • ). When θ F is smaller than 45 • , the projected contour of the top view of the bubble is mainly affected by the length of axis B, so axis A is rotated to the vertical direction; on the other hand, when θ F is greater than 45 • , the projected contour of the top view of the bubble is mostly affected by axis A, so axis A is rotated to the horizontal direction. The rotation formula is given by: where, (x F , z F ) and (x F , z F ) represent the coordinates of the bubble projection contour before rotation and after rotation, respectively; θ F is the angle of rotation. Then, the magnification factor of the front view is calculated by: where, γ F is the magnification factor of the front view; S F and s F represent the maximum width of the front image and that of the corrected front image, respectively. Meanwhile, the coincident steps were repeated to correct the projected contour of the side image, and the magnification factor of the side view (γ S ) was then figured. Finally, in order to make the bubble contours conform to physical and geometric relationship, the top projection contour needed to be corrected using the magnification of the front view and side view as follows: where, (x T , y T ) are the coordinates of the uncorrected top view image, and (x T , y T ) denote the coordinates of the corrected top view image. After the top view image was corrected, the angle between the long axis (or the short axis) and the horizontal direction was extracted and that was regarded as the azimuth angle (β T ). where, γF is the magnification factor of the front view; SF and sF represent the maximum width of the front image and that of the corrected front image, respectively. Meanwhile, the coincident steps were repeated to correct the projected contour of the side image, and the magnification factor of the side view (γS) was then figured. Finally, in order to make the bubble contours conform to physical and geometric relationship, the top projection contour needed to be corrected using the magnification of the front view and side view as follows: where, (xT, yT) are the coordinates of the uncorrected top view image, and ( T x , T y ) denote the coordinates of the corrected top view image. After the top view image was corrected, the angle between the long axis (or the short axis) and the horizontal direction was extracted and that was regarded as the azimuth angle (βT).

Calculating Parameters of Horizontal Ellipse Slices
Since the dimension of the single bubble is much smaller than the distance between the camera and the bubble, it can be reasonable to assume that the incident light is parallel light, so the bubble is evenly divided into n pieces of ellipse slices along the height direction. Further, both of the azimuth angle (βT) and the projection parameter, α (αF or αS), of the ellipse slices were assumed to be constant at each layer. The major and minor semi-axes of the ith layer of the ellipse at a height of Hi were indicated by ai and bi, respectively, and the projection widths in the front and side views were

Calculating Parameters of Horizontal Ellipse Slices
Since the dimension of the single bubble is much smaller than the distance between the camera and the bubble, it can be reasonable to assume that the incident light is parallel light, so the bubble is evenly divided into n pieces of ellipse slices along the height direction. Further, both of the azimuth angle (β T ) and the projection parameter, α (α F or α S ), of the ellipse slices were assumed to be constant at each layer. The major and minor semi-axes of the ith layer of the ellipse at a height of H i were indicated by a i and bi, respectively, and the projection widths in the front and side views were denoted by l Fi and l Si , respectively ( Figure 6). The equation of the ith layer of the ellipse projected onto the top view can be written by: where, A i = Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 18 denoted by lFi and lSi, respectively ( Figure 6). The equation of the ith layer of the ellipse projected onto the top view can be written by: where, 22 22 cos sin Two tangent points are derived from the two tangent lines parallel with the x-axis on layer i, and the chord length of the line passing through the two tangent points and the ellipse is equal to the width of projection onto the side view (lSi) as shown in Figure 6. The formula for lSi is expressed as: where, S  is a projection parameter associated with the angle between the chord length and horizontal line. By rearranging, these parameters can be related as follows: The tangent coordinates of the tangent line, the slope of which is zero since it is parallel with the x-axis, are defined by (x, x tan αS), where: Equation (11) could be simplified as follows: where, cos sin tan sin cos sin tan cos By combining Equations (10) and (12), one can obtain: Two tangent points are derived from the two tangent lines parallel with the x-axis on layer i, and the chord length of the line passing through the two tangent points and the ellipse is equal to the width of projection onto the side view (l Si ) as shown in Figure 6. The formula for l Si is expressed as: where, α S is a projection parameter associated with the angle between the chord length and horizontal line. By rearranging, these parameters can be related as follows: By rearranging Equation (9) into a dimensionless equation, where a 0 = a i l Si The tangent coordinates of the tangent line, the slope of which is zero since it is parallel with the x-axis, are defined by (x, x tan α S ), where: Equation (11) could be simplified as follows: Appl. Sci. 2020, 10, 5803 9 of 18 By combining Equations (10) and (12), one can obtain: Equation (13) shows that the dimensionless parameters a 0 and b 0 are only associated with the azimuth angle (β T ), so the projection parameter (α S ) of an arbitrary ellipse slice is constant. Both the long semi-axis (a i ) and short semi-axis (b i ) of the ellipse slice in the ith layer can be calculated by the following formula: Using the same procedure, similar results can also be obtained by using the projection width of the front view, and the corresponding dimensionless parameters of the ellipse can be calculated as: The long semi-axis (a i ) and short semi-axis (b i ) in layer i are defined using the following equation: All the ellipse slices were then stacked to form the 3D shape of the single bubble after traversing each layer of the ellipse slice. It is worth noting that when β T approached 90 • inseparably, the solutions to Equations (13) and (15) were not practical in the reconstruction process under these conditions; therefore, it is reasonable to assume that both axes of each ellipse slice are equal to the corresponding projection width to avoid the odd points of Equations (13) and (15).

Bubble Samples
Many studies used bubble-like objects to test the accuracy of reconstruction [18,23]. Therefore, three different types of bubble-like object (hereinafter referred to as bubble) were produced using 3D printing technology as displayed in Figure 7, to evaluate the measurement accuracy of several reconstruction methods discussed above. The first type was a sphere, where all the semi-axes have equal lengths, i.e., a = b = c; this type is generally used in ideal situations [31,32]. The second type was an ellipsoid, in which at least two of the three semi-axes (a, b, and c) were unequal; this type is also suitable for describing the bubble in a shear flow [33,34]. The third type resembled an inverted teardrop, which is mostly developed in viscose or shear thinning fluids [1,4]. Since the finite acrylic box existed wall influence, which led to uncertainties and different a, b, and c values, a series of objects with identical a and c but different b values were generated. An inverted droplet was produced using equal values of a and b, denoting a circular horizontal cross-section. Furthermore, the accuracy of 3D printing was ±0.01 mm. Hence, the total measurement accuracy calculated by Equation (1) was ±0.46 mm. The accuracy of volume, tested by the densimeter AU-300PM, was ±1 mm 3 . The locations of cameras were stationary during experiments. Thus, the bubble parameters can be measured (Table 3) and the objects are shown in Figure 8.

Effects of Posture Correction
The bubble volume was used to measure the accuracy of the reconstruction results. The difference between the calculated volume and the actual volume in the form of mean relative error (MRE) was calculated during experiments. The standard deviation of the calculated results was represented by error bars. All the examples were tested three times and the rotation angles around x-axis, y-axis, and z-axis of each test case were generated randomly. Figure 9 compares the accuracy of the reconstruction methods (in the form of MRE) with/without posture correction, where the sphere (sample No. 12), ellipsoid (sample No. 7), and inverted teardrop (sample No. 13) are analyzed using the methods in this paper, Fujiwara et al. [21], and Bian et al. [20]. Since the spheres have an equivalent length in whichever direction, posture correction does not theoretically influence the reconstruction, leading to slight distinction between

Effects of Posture Correction
The bubble volume was used to measure the accuracy of the reconstruction results. The difference between the calculated volume and the actual volume in the form of mean relative error (MRE) was calculated during experiments. The standard deviation of the calculated results was represented by error bars. All the examples were tested three times and the rotation angles around x-axis, y-axis, and z-axis of each test case were generated randomly. Figure 9 compares the accuracy of the reconstruction methods (in the form of MRE) with/without posture correction, where the sphere (sample No. 12), ellipsoid (sample No. 7), and inverted teardrop (sample No. 13) are analyzed using the methods in this paper, Fujiwara et al. [21], and Bian et al. [20]. Since the spheres have an equivalent length in whichever direction, posture correction does not theoretically influence the reconstruction, leading to slight distinction between the corrected and uncorrected data. On the other hand, the ellipsoids reflect the degree of deformation in the horizontal direction, while the inverted teardrops show the degree of deformation in the vertical direction. Therefore, by comparing the data on the ellipsoids and the inverted teardrops, it can logically be found that posture correction significantly enhances the reconstruction accuracy of the method in this paper, but the effects on methods in Fujiwara et al. [21] and Bian et al. [20] are slight. This indicates that the posture correction is undoubtedly required in our method for the reconstruction of bubbles having serious deformation.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 18 the corrected and uncorrected data. On the other hand, the ellipsoids reflect the degree of deformation in the horizontal direction, while the inverted teardrops show the degree of deformation in the vertical direction. Therefore, by comparing the data on the ellipsoids and the inverted teardrops, it can logically be found that posture correction significantly enhances the reconstruction accuracy of the method in this paper, but the effects on methods in Fujiwara et al. [21] and Bian et al. [20] are slight. This indicates that the posture correction is undoubtedly required in our method for the reconstruction of bubbles having serious deformation.

Effects of Different Shapes
It can be seen in Figure 10 that the MRE calculated for the three methods is approximately less than 3% when the bubble resembles a sphere, which indicates that the three methods can well reconstruct spherical bubbles. Nevertheless, when the bubble is deformed in the horizontal direction, such as ellipsoidal bubbles, the accuracy of the reconstruction methods changes progressively in the order this paper > Fujiwara et al. [21] > Bian et al. [20], that is, the modified stacked ellipse method leads to the smallest mean relative error in reconstructing the bubble. In fact, since the method in Fujiwara et al. [21] utilizes the projection width onto the front and side views as the long and short axes of the ellipse on each ellipse slice respectively, when the major axis of the ellipse is not equal to the projection width, that is, βT is deviated from either 0 or 90°, the method in Fujiwara et al. [21] leads to a larger error in reconstructing the bubble (this phenomenon could also been seen in Section 4.5). Furthermore, the assumption of Bian et al. [20] that a bubble has a minimum rotation angle around z-axis was overturned because the azimuth angle βT always exists, and the bubble's highest point is not always located at the z-axis. Therefore, Bian et al. [20] overestimate the bubble, especially when the bubble is not elliptical. In this way, the method in Bian et al. [20] leads to a greater error in reconstructing the bubble compared with the method in Fujiwara et al. [21]. When the deformation of the bubble is increased in the vertical direction, that is, the bubble has an inverted-teardrop shape. The method in this paper shows a surprising degree of accuracy, and the method in Fujiwara et al. [21] is also far superior to the method in Bian et al. [20], which is attributed to the fact that the methods in this paper and Fujiwara et al. [21] stack ellipse slices along the vertical direction, causing a minimal impact on the reconstruction accuracy when vertical deformation increases. The method in Bian et al. [20] assumes that the vertical projection is oval, which leads to an extremely large MRE in reconstructing an inverted teardrop-shaped bubble.

Effects of Different Shapes
It can be seen in Figure 10 that the MRE calculated for the three methods is approximately less than 3% when the bubble resembles a sphere, which indicates that the three methods can well reconstruct spherical bubbles. Nevertheless, when the bubble is deformed in the horizontal direction, such as ellipsoidal bubbles, the accuracy of the reconstruction methods changes progressively in the order this paper > Fujiwara et al. [21] > Bian et al. [20], that is, the modified stacked ellipse method leads to the smallest mean relative error in reconstructing the bubble. In fact, since the method in Fujiwara et al. [21] utilizes the projection width onto the front and side views as the long and short axes of the ellipse on each ellipse slice respectively, when the major axis of the ellipse is not equal to the projection width, that is, β T is deviated from either 0 or 90 • , the method in Fujiwara et al. [21] leads to a larger error in reconstructing the bubble (this phenomenon could also been seen in Section 4.5). Furthermore, the assumption of Bian et al. [20] that a bubble has a minimum rotation angle around z-axis was overturned because the azimuth angle β T always exists, and the bubble's highest point is not always located at the z-axis. Therefore, Bian et al. [20] overestimate the bubble, especially when the bubble is not elliptical. In this way, the method in Bian et al. [20] leads to a greater error in reconstructing the bubble compared with the method in Fujiwara et al. [21]. When the deformation of the bubble is increased in the vertical direction, that is, the bubble has an inverted-teardrop shape. The method in this paper shows a surprising degree of accuracy, and the method in Fujiwara et al. [21] is also far superior to the method in Bian et al. [20], which is attributed to the fact that the methods in this paper and Fujiwara et al. [21] stack ellipse slices along the vertical direction, causing a minimal impact on the reconstruction accuracy when vertical deformation increases. The method in Bian et al. [20] assumes that the vertical projection is oval, which leads to an extremely large MRE in reconstructing an inverted teardrop-shaped bubble.

Effects of Deformation
The above analysis revealed that the degree of bubble deformation has a significant effect on the reconstruction accuracy. Therefore, the following ellipsoids (samples No. [1][2][3][4][5][6][7][8][9][10][11] were used as examples to analyze the influences of the horizontal deformation of the bubble, in terms of the ratio of the vertical axis to the horizontal axis, i.e., (b/a), and the sphericity of the bubble on the accuracy of the reconstruction algorithms. Figure 11 shows that the MRE values of the three methods are similar when the ratio of the vertical axis to the horizontal axis, (b/a), is greater than 0.5 and confirms that the method in this paper performs better than the two other methods with a reconstruction error smaller than 6%. Reducing the ratio of the vertical axis to the horizontal axis in the range of 0 to 0.5 increases the mean relative error obtained by the three methods. Moreover, the mean relative error in reconstructing the bubble by the method in this paper is not more than 24%, which denotes a satisfying degree of accuracy. Nonetheless, the accuracy of the method in Fujiwara et al. [21] is closely similar to that of the method in Bian et al. [20] since all the analyzed samples are standard ellipsoids, the vertical projection of which is oval. Therefore, the capability of the method in Fujiwara et al. [21] to resist vertical deformation does not give it an advantage over the method in Bian et al. [20] in this case.
where, SS denotes the surface area of a sphere having the same volume as the object, and SO stands for the surface area of the object. Figure 12 demonstrates that the accuracy of the three methods in calculating the bubble volume, or rather reconstructing the bubble, is not much different when the sphericity is greater than 0.90; however, the method in this paper is slightly more accurate than the two other methods.

Effects of Deformation
The above analysis revealed that the degree of bubble deformation has a significant effect on the reconstruction accuracy. Therefore, the following ellipsoids (samples No. [1][2][3][4][5][6][7][8][9][10][11] were used as examples to analyze the influences of the horizontal deformation of the bubble, in terms of the ratio of the vertical axis to the horizontal axis, i.e., (b/a), and the sphericity of the bubble on the accuracy of the reconstruction algorithms. Figure 11 shows that the MRE values of the three methods are similar when the ratio of the vertical axis to the horizontal axis, (b/a), is greater than 0.5 and confirms that the method in this paper performs better than the two other methods with a reconstruction error smaller than 6%. Reducing the ratio of the vertical axis to the horizontal axis in the range of 0 to 0.5 increases the mean relative error obtained by the three methods. Moreover, the mean relative error in reconstructing the bubble by the method in this paper is not more than 24%, which denotes a satisfying degree of accuracy. Nonetheless, the accuracy of the method in Fujiwara et al. [21] is closely similar to that of the method in Bian et al. [20] since all the analyzed samples are standard ellipsoids, the vertical projection of which is oval. Therefore, the capability of the method in Fujiwara et al. [21] to resist vertical deformation does not give it an advantage over the method in Bian et al. [20] in this case.

Effects of Deformation
The above analysis revealed that the degree of bubble deformation has a significant effect on the reconstruction accuracy. Therefore, the following ellipsoids (samples No. [1][2][3][4][5][6][7][8][9][10][11] were used as examples to analyze the influences of the horizontal deformation of the bubble, in terms of the ratio of the vertical axis to the horizontal axis, i.e., (b/a), and the sphericity of the bubble on the accuracy of the reconstruction algorithms. Figure 11 shows that the MRE values of the three methods are similar when the ratio of the vertical axis to the horizontal axis, (b/a), is greater than 0.5 and confirms that the method in this paper performs better than the two other methods with a reconstruction error smaller than 6%. Reducing the ratio of the vertical axis to the horizontal axis in the range of 0 to 0.5 increases the mean relative error obtained by the three methods. Moreover, the mean relative error in reconstructing the bubble by the method in this paper is not more than 24%, which denotes a satisfying degree of accuracy. Nonetheless, the accuracy of the method in Fujiwara et al. [21] is closely similar to that of the method in Bian et al. [20] since all the analyzed samples are standard ellipsoids, the vertical projection of which is oval. Therefore, the capability of the method in Fujiwara et al. [21] to resist vertical deformation does not give it an advantage over the method in Bian et al. [20] in this case. Sphericity (s), as a parameter measuring the degree of deviation of an object from a sphere, is defined by: where, SS denotes the surface area of a sphere having the same volume as the object, and SO stands for the surface area of the object. Figure 12 demonstrates that the accuracy of the three methods in calculating the bubble volume, or rather reconstructing the bubble, is not much different when the sphericity is greater than 0.90; however, the method in this paper is slightly more accurate than the two other methods. Sphericity (s), as a parameter measuring the degree of deviation of an object from a sphere, is defined by: where, S S denotes the surface area of a sphere having the same volume as the object, and S O stands for the surface area of the object. Figure 12 demonstrates that the accuracy of the three methods in calculating the bubble volume, or rather reconstructing the bubble, is not much different when the sphericity is greater than 0.90; however, the method in this paper is slightly more accurate than the two other methods.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 18 Therefore, based on the above discussion, the method in this paper is more accurate and has better applicability when the bubble is horizontally or vertically deformed.

Effects of Azimuth Angle
The effect of azimuth angle on the accuracy of the reconstruction methods is studied by reconstructing three ellipsoids (samples No. 1, 3, and 7), as delineated in Figure 13. The mean relative error in the bubble volume calculated by the three methods, which denotes the accuracy of the methods in reconstructing the bubble, is similar when the azimuth angle is close to either 0 or 90°. However, when the azimuth angle gradually deviates from either 0 or 90°, the mean relative errors calculated by both the method in Fujiwara et al. [21] and Bian et al. [20] are obviously larger than the MRE obtained by the method in this paper. Figure 13 shows the methods of Fujiwara et al. [21] and Bian et al. [20] were more sensitive to bubble shape and rotation angles. Besides, the method in this paper was more stable in measuring bubble volume. It is notable that the bubble shapes and rotation angles had little impact on the volume measurement accuracy using the method in this paper. The ellipsoid (No.1) in 3rd graph of Figure 13 was visualized in Figure 14 as an example, with their corresponding projection images for front view, side view, and top view. As shown in case 1 to 3 of Figure 14, the projection contours always overstated the bubble size. If bubble rotated around z axis, the front view and side view projection contours would become much larger. The reconstruction methods based on Fujiwara et al. [21] and Bian et al. [20], by this token, would have larger uncertainties and unstable predictions. This scenario got significant when the azimuth angle was approximate to 45°. Therefore, based on the above discussion, the method in this paper is more accurate and has better applicability when the bubble is horizontally or vertically deformed.

Effects of Azimuth Angle
The effect of azimuth angle on the accuracy of the reconstruction methods is studied by reconstructing three ellipsoids (samples No. 1, 3, and 7), as delineated in Figure 13. The mean relative error in the bubble volume calculated by the three methods, which denotes the accuracy of the methods in reconstructing the bubble, is similar when the azimuth angle is close to either 0 or 90 • . However, when the azimuth angle gradually deviates from either 0 or 90 • , the mean relative errors calculated by both the method in Fujiwara et al. [21] and Bian et al. [20] are obviously larger than the MRE obtained by the method in this paper. Figure 13 shows the methods of Fujiwara et al. [21] and Bian et al. [20] were more sensitive to bubble shape and rotation angles. Besides, the method in this paper was more stable in measuring bubble volume. It is notable that the bubble shapes and rotation angles had little impact on the volume measurement accuracy using the method in this paper.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 18 Therefore, based on the above discussion, the method in this paper is more accurate and has better applicability when the bubble is horizontally or vertically deformed.

Effects of Azimuth Angle
The effect of azimuth angle on the accuracy of the reconstruction methods is studied by reconstructing three ellipsoids (samples No. 1, 3, and 7), as delineated in Figure 13. The mean relative error in the bubble volume calculated by the three methods, which denotes the accuracy of the methods in reconstructing the bubble, is similar when the azimuth angle is close to either 0 or 90°. However, when the azimuth angle gradually deviates from either 0 or 90°, the mean relative errors calculated by both the method in Fujiwara et al. [21] and Bian et al. [20] are obviously larger than the MRE obtained by the method in this paper. Figure 13 shows the methods of Fujiwara et al. [21] and Bian et al. [20] were more sensitive to bubble shape and rotation angles. Besides, the method in this paper was more stable in measuring bubble volume. It is notable that the bubble shapes and rotation angles had little impact on the volume measurement accuracy using the method in this paper. The ellipsoid (No.1) in 3rd graph of Figure 13 was visualized in Figure 14 as an example, with their corresponding projection images for front view, side view, and top view. As shown in case 1 to 3 of Figure 14, the projection contours always overstated the bubble size. If bubble rotated around z axis, the front view and side view projection contours would become much larger. The reconstruction methods based on Fujiwara et al. [21] and Bian et al. [20], by this token, would have larger uncertainties and unstable predictions. This scenario got significant when the azimuth angle was approximate to 45°. The ellipsoid (No.1) in 3rd graph of Figure 13 was visualized in Figure 14 as an example, with their corresponding projection images for front view, side view, and top view. As shown in case 1 to 3 of Figure 14, the projection contours always overstated the bubble size. If bubble rotated around z axis, the front view and side view projection contours would become much larger. The reconstruction methods based on Fujiwara et al. [21] and Bian et al. [20], by this token, would have larger uncertainties and unstable predictions. This scenario got significant when the azimuth angle was approximate to 45 • .

Verification by Gas Injection Test
To further verify the accuracy of the method in this paper, a gas injection test was conducted using the experimental apparatus in Section 2. During this experiment, the bubble images were taken by three cameras, the gas flow rate was set as 1 sccm, and the gas pressure was recorded by pressure sensor sequentially. As we can see in Figure 15, the pressure increased steadily until it met a peak, and then dropped ephemerally. After that, it maintains a stable horizontal level. The bubble appeared soon after the pressure reached peak and it grew during the period of pressure drop and the stable pressure level. Images at four different time nodes (A-D) were stitched together to visualize the bubble growth. The experimental bubble volumes were calculated by following steps based on Van der Waals (1873). The initial condition (point O0) could be written as: where, P0 is the initial pipe pressure, and the n0 is initial mole numbers inside the pipe, Vpipe is the pipe volume.
By combining Equations (18) and (19), both P0 and n0 were solved. For the period of pressure drop (point O1 to O2) in Figure 15:

Verification by Gas Injection Test
To further verify the accuracy of the method in this paper, a gas injection test was conducted using the experimental apparatus in Section 2. During this experiment, the bubble images were taken by three cameras, the gas flow rate was set as 1 sccm, and the gas pressure was recorded by pressure sensor sequentially. As we can see in Figure 15, the pressure increased steadily until it met a peak, and then dropped ephemerally. After that, it maintains a stable horizontal level. The bubble appeared soon after the pressure reached peak and it grew during the period of pressure drop and the stable pressure level. Images at four different time nodes (A-D) were stitched together to visualize the bubble growth.

Verification by Gas Injection Test
To further verify the accuracy of the method in this paper, a gas injection test was conducted using the experimental apparatus in Section 2. During this experiment, the bubble images were taken by three cameras, the gas flow rate was set as 1 sccm, and the gas pressure was recorded by pressure sensor sequentially. As we can see in Figure 15, the pressure increased steadily until it met a peak, and then dropped ephemerally. After that, it maintains a stable horizontal level. The bubble appeared soon after the pressure reached peak and it grew during the period of pressure drop and the stable pressure level. Images at four different time nodes (A-D) were stitched together to visualize the bubble growth. The experimental bubble volumes were calculated by following steps based on Van der Waals (1873). The initial condition (point O0) could be written as: where, P0 is the initial pipe pressure, and the n0 is initial mole numbers inside the pipe, Vpipe is the pipe volume.
By combining Equations (18) and (19), both P0 and n0 were solved. For the period of pressure drop (point O1 to O2) in Figure 15: The experimental bubble volumes were calculated by following steps based on Van der Waals (1873). The initial condition (point O 0 ) could be written as: where, P 0 is the initial pipe pressure, and the n 0 is initial mole numbers inside the pipe, V pipe is the pipe volume. The Van der Waals parameters a = 0.1408 Pa · m 6 · mol −2 , b = 3.912 × 10 −5 m 5 · mol −1 , and the ideal gas constant R = 8.31441 ± 0.00026 J/(mol · K). The thermodynamic temperature K inside the laboratory was 288 K. For a moment (point O') where the bubble had not generated: By combining Equations (18) and (19), both P 0 and n 0 were solved. For the period of pressure drop (point O 1 to O 2 ) in Figure 15: where, P 1 is the beginning pressure of pressure drop, and P 2 is the ending pressure of pressure drop, V b2 is the calculated bubble volume at the point O 2 , n 2 is the mole numbers at the point O 2 : Then, the bubble volume (V b ) after O 2 are: where, t x is the experimental time. A series of experiment pictures (corresponding to points A-D respectively in Figure 15) were taken based on the setup, some of which were shown as example in Figure 16. To further testify the accuracy of the method proposed above, they were reconstructed in Figure 17. It can be seen that the bubbles developed vertically and the projections for horizontal direction were analogue. The bubble shape was close to cone at the beginning, and then developed to ellipsoid quickly. The tip of the cusp emerged as bubble volume increased, and eventually evolved into shape of inverted teardrop.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 18 where, P1 is the beginning pressure of pressure drop, and P2 is the ending pressure of pressure drop, Vb2 is the calculated bubble volume at the point O2, n2 is the mole numbers at the point O2: Then, the bubble volume (Vb) after O2 are: where, tx is the experimental time. A series of experiment pictures (corresponding to points A-D respectively in Figure 15) were taken based on the setup, some of which were shown as example in Figure 16. To further testify the accuracy of the method proposed above, they were reconstructed in Figure 17. It can be seen that the bubbles developed vertically and the projections for horizontal direction were analogue. The bubble shape was close to cone at the beginning, and then developed to ellipsoid quickly. The tip of the cusp emerged as bubble volume increased, and eventually evolved into shape of inverted teardrop.  The bubble volume measured by the methods in this paper, in Fujiwara et al. [21] and in Bian et al. [20], are shown in Figure 18. The construction results of the method in this paper fits well with the theoretical results and the MRE of bubble volume later O2, which was 1.79%. Similarly, the where, P1 is the beginning pressure of pressure drop, and P2 is the ending pressure of pressure drop, Vb2 is the calculated bubble volume at the point O2, n2 is the mole numbers at the point O2: Then, the bubble volume (Vb) after O2 are: where, tx is the experimental time. A series of experiment pictures (corresponding to points A-D respectively in Figure 15) were taken based on the setup, some of which were shown as example in Figure 16. To further testify the accuracy of the method proposed above, they were reconstructed in Figure 17. It can be seen that the bubbles developed vertically and the projections for horizontal direction were analogue. The bubble shape was close to cone at the beginning, and then developed to ellipsoid quickly. The tip of the cusp emerged as bubble volume increased, and eventually evolved into shape of inverted teardrop.  The bubble volume measured by the methods in this paper, in Fujiwara et al. [21] and in Bian et al. [20], are shown in Figure 18. The construction results of the method in this paper fits well with the theoretical results and the MRE of bubble volume later O2, which was 1.79%. Similarly, the The bubble volume measured by the methods in this paper, in Fujiwara et al. [21] and in Bian et al. [20], are shown in Figure 18. The construction results of the method in this paper fits well with the theoretical results and the MRE of bubble volume later O 2 , which was 1.79%. Similarly, the construction results of the method in Fujiwara et al. [21] also fits the theoretical results perfectly, and the MRE of bubble volume later O 2 was 1.74%. The construction results of the method in Bian et al. [20] fits the theoretical results at the beginning; however, the error arose as time went by. In the method of Bian et al. [20], the MRE of bubble volume was 21.2%. As observed in gas injection test, the horizontal deformation of bubble was negligible due to the same effect of wall effect. Therefore, the construction results in this paper were almost the same as that in Fujiwara et al. [21]. Furthermore" the vertical shape of bubble developed from ellipses to inverted-droplet, causing the construction accuracy high firstly and reduced later. It indicated that the method in this paper per have credible consistency with experimental results where bubbles had larger vertical deformation and little horizontal deformation.

Conclusions
In the present work, we developed a precise method based on orthographic images to reconstruct a 3D bubble shape of a single bubble, which is elliptical in the horizontal direction and has an arbitrary shape in the vertical direction. The proposed reconstruction method demonstrated a surprising degree of accuracy compared with the methods in Fujiwara et al. [21] and Bian et al. [20], and its superiority was notable when the bubble deformation was increased.
The method in this paper differed from Fujiwara et al. [21] in two points. Firstly, the posture correction reduced the tilt problem caused by bubble oscillation. Second, the projection relationship and the geometric information were used to calculate the major and minor semi-axes of each horizontal ellipse slices. Therefore, instead of using the projection width directly, the method in this paper reduced the negative impact of horizontal deformation and azimuth angle on the reconstruction accuracy.
According to the above discussion, the following conclusions can be drawn: (1) The posture correction helps to improve the reconstruction accuracy in this paper, especially for the bubble with a horizontal cross-section remarkably deviating from a circle. (2) The dimensionless parameters of an ellipse in horizontal section (major and minor semi-axes) are independent of the projection width. By using the dimensionless elliptic parameters, the actual major and minor semi-axes of the horizontal elliptical section of each layer can be quickly calculated. (3) The method in Bian et al. [20] is more suitable for reconstructing a bubble similar to an ellipsoid because both of the horizontal and vertical deformations inversely impact its accuracy. The vertical deformation has a little influence on the method in Fujiwara et al. [21], but the horizontal deformation contributes significantly to its accuracy, especially when the ratio of the vertical axis to the horizontal axis is smaller than 0.5 or when the azimuth angle deviates greatly from either 0 or 90°. (4) The method in this paper successfully improves the flaws in methods in Bian et al. [20] and Fujiwara et al. [21]. The method in this paper can better reconstruct the bubble with severe horizontal deformation; in addition, similarly to the method in Fujiwara et al. [21], the vertical deformation has a negligible influence on the reconstruction accuracy of the method in this paper, leading to its higher practicality.

Conclusions
In the present work, we developed a precise method based on orthographic images to reconstruct a 3D bubble shape of a single bubble, which is elliptical in the horizontal direction and has an arbitrary shape in the vertical direction. The proposed reconstruction method demonstrated a surprising degree of accuracy compared with the methods in Fujiwara et al. [21] and Bian et al. [20], and its superiority was notable when the bubble deformation was increased.
The method in this paper differed from Fujiwara et al. [21] in two points. Firstly, the posture correction reduced the tilt problem caused by bubble oscillation. Second, the projection relationship and the geometric information were used to calculate the major and minor semi-axes of each horizontal ellipse slices. Therefore, instead of using the projection width directly, the method in this paper reduced the negative impact of horizontal deformation and azimuth angle on the reconstruction accuracy.
According to the above discussion, the following conclusions can be drawn: (1) The posture correction helps to improve the reconstruction accuracy in this paper, especially for the bubble with a horizontal cross-section remarkably deviating from a circle. (2) The dimensionless parameters of an ellipse in horizontal section (major and minor semi-axes) are independent of the projection width. By using the dimensionless elliptic parameters, the actual major and minor semi-axes of the horizontal elliptical section of each layer can be quickly calculated. (3) The method in Bian et al. [20] is more suitable for reconstructing a bubble similar to an ellipsoid because both of the horizontal and vertical deformations inversely impact its accuracy. The vertical deformation has a little influence on the method in Fujiwara et al. [21], but the horizontal deformation contributes significantly to its accuracy, especially when the ratio of the vertical axis to the horizontal axis is smaller than 0.5 or when the azimuth angle deviates greatly from either 0 or 90 • . (4) The method in this paper successfully improves the flaws in methods in Bian et al. [20] and Fujiwara et al. [21]. The method in this paper can better reconstruct the bubble with severe horizontal deformation; in addition, similarly to the method in Fujiwara et al. [21], the vertical deformation has a negligible influence on the reconstruction accuracy of the method in this paper, leading to its higher practicality. (5) The two methods in Bian et al. [20] and Fujiwara et al. [21] are extremely sensitive to the orientation of the bubble since the azimuth angle of the horizontal elliptical slice strongly affects the reconstruction accuracy; however, it does not impact the method in this paper.
The bubble-shaped samples used in our experiments were considered to be regular bubbles since the bubbles generally representing a larger volume, less distortion, and a regular shape in Newtonian or non-Newtonian fluids with higher viscosity are similar to regular bubbles. Therefore, the reconstruction method developed herein is undoubtedly suitable for such a kind of bubbles. However, the application of the proposed method to reconstructing bubbles with a large horizontal and/or vertical distortion must be further validated. In this way, further exploration needs to be done to study the accuracy of this method when bubbles have both vertical and horizontal deformation. The algorithm should be sufficiently replenished on the problem that cameras are not ideal when perpendicular to each other.

Conflicts of Interest:
The authors declare no conflict of interest.