3D Measurements of a Two-Phase Flow Inside an Optical Cylinder Based on Full-Field Cross-Interface Computed Tomography

Three-dimensional (3D) tomographic reconstruction in confined-space requires a mapping relationship which considers the refraction distortion caused by optical walls. In this work, a tomography method, namely full-field cross-interface computed tomography (FCICT), is proposed to solve confine-space problems. The FCICT method utilizes Snell’s law and reverse ray-tracing to analytically correct imaging distortion and establishes the mapping relationship from 3D measurement domain to 2D images. Numerical phantom study is first employed to validate the FCICT method. Afterwards, the FCICT is applied on the experimental reconstruction of an illuminated two-phase jet flow which is initially generated inside an optical cylinder and then gradually moves outside. The comparison between accurately reconstructed liquid jet by FCICT and coarse result by traditional open space tomography algorithm provides a practical validation of FCICT. Based on the 3D liquid jet reconstructions at different time sequences, the distributions of surface velocity and 3D curvatures are calculated, and their correspondences are systematically analyzed. It is found that the velocity of a surface point is positively correlated with the mean curvature at the same point, which indicates the concavity/convexity of liquid jet surface is possibly in accordance with the surface velocity. Moreover, the surface velocity presents monotonical increasing trend with larger Gaussian curvature for elliptic surface points only, due to the dominated Brownian motion as the liquid jet develops.


Introduction
3D tomography has become an irreplaceable approach in flow and combustion researches (Yang et al. 2021;Hwang et al. 2020). The basic principle of 3D tomography is to reconstruct the 3D optical field by mathematical inversion algorithm using a set of 2D images of the 3D optical field and the mapping relationship between them. To reconstruct the 3D distribution of different practical optical signals, the tomography problem can be divided into open-space and confined-space problems. For open-space tomography problems, the signal from optical field is not blocked, refracted or reflected before it is captured by the optical detectors. Hence, the projecting relationships between signal field and detectors can be determined easily by pin-hole model or other imaging models. Past researches have made solid progress on open-space tomography problems. For example, the topology and surface properties of a matrix burner are analyzed according to 3D flame reconstruction from 3D tomography method (Floyd and Kempf 2011). Volumetric laser-induced fluorescence (VLIF) based on tomography was developed to measure the 4D structures of dynamic flames and temperature fields in turbulent flames (Ma et al. 2017;Halls et al. 2018). Traditional 2D particle image velocimetry (PIV) and 3D tomography were combined to measure the 3D velocity distribution in flow fields (Ebi and Clemens 2016;Liu et al. 2018). The triangle-mesh-based method and tomography were combined to calculate the 3D surface curvatures (i.e., mean curvature and Gaussian curvatures) of turbulent swirl flames (Yu et al. 2020). For confined-space tomography problems, the optical signal field is enclosed by the facility walls while the imaging system is set outside. Some portions of the facility walls are made of transparent materials to allow the access of optical signal, such as optical internal combustion engines (Zhang et al. 2019) and wind/water tunnels (Cao et al. 2021;Zhang and Vanierschot 2021). Examples of 3D tomography approaches in these facilities include the burning velocity, topology and winked structure measurements of turbulent flames in a 1 MPa combustion chamber (Zhang et al. 2019), and the flame ignition dynamics measurements in a Mach 2 combustor (Ma et al. 2015). It is worth noting that the optical window surfaces are flat in abovementioned examples. Under such situation, the imaging distortion appears to be integral translation for all locations on the projection, which can be easily fixed by integrally shifting the 2D signal distribution in each projection during the calibration process (Ma et al. 2016a).
When the facility has curved walls, such as the transparent cylinder in optical engines, the distortion of projection is nonlinear and cannot be fixed only by integral shifting. To solve such problems, past efforts primarily focused on developing numerical corrections on the imaging process across refractive medium. For example, empirical polynomial models were developed to solve the distortion problems (Prasad 2000;Soloff et al. 1997;Falkhytten 2018). The main drawback of polynomial models is that they are generated based on specific optical wall parameters (e.g., curvatures, refractive index, etc.), thus cannot be utilized universally. Some studies calculated the light path after double-refraction on the optical wall by ray-tracing method based on pinhole model (Kotowski 1988;Mulsow 2010). However, these studies employed a large number of iterations to accurately determine the light path after refraction, which is time-consuming. To avoid iterations, reversed ray-tracing method was developed to obtain the mapping relationship from camera chips to 3D optical signal field (Mulsow 2010;Belden 2013;Liu et al. 2019). However, these approaches cannot solve the forward ray-tracing relationship from discretized 3D target (i.e., voxels) in the measurement domain to the pixels on the sensors. Based on the understanding of past efforts, we proposed an improved tomography method, namely full-field cross-interface computed tomography (FCICT) (Ling et al. 2021(Ling et al. , 2020, to solve the confine-space problems precisely and efficiently. The reverse ray-tracing method and Snell's refraction law are used in the FCICT to establish the distorted mapping relationship. The correctness and necessity of FCICT are validated by both simulation and controlled laminar-flow cone-flame experiments. Based on the algorithm development and controlled experiment validation of the FCICT, the method is applied to perform the 3D measurements of two-phase flow (a gaseous jet seeded with micro-sized water droplets) inside an optical cylinder. The two-phase flow is illuminated by a continuous wave laser, and a high-speed camera is employed to the evolution of the flow patterns. Based on the FCICT reconstruction sequences, the 3D curvature and diffusion velocity distribution on the flow surface are calculated and analyzed, as a representative results of 3D tomography. In the rest of this paper, Sect. 2 describes the details of the FCICT algorithm with numerical validation based on simulated 3D phantom. Section 3 describes the experimental setup of the two-phase flow 3D measurements, and then compares the flow reconstructions obtained by FCICT and traditional open-space tomography method to prove the accuracy of FCICT. Section 4 demonstrates and analyzes the 3D curvature and velocity distribution achieved from reconstructed flow sequences. Finally, Sect. 5 concludes the paper.

The FCICT Algorithm Description and Validation
This section describes the algorithm details of FCICT algorithm used to solve the practical confined-space tomographic problem. The confined-space is accomplished by installing the cylinder of an optical internal combustion engine around the measurement domain. Similar to traditional tomography algorithms, the FCICT algorithm consists of two parts: (1) the establishment of the point spread function (PSF) matrix that denotes the signal mapping relationship between voxels in the measurement domain and pixels on the sensor arrays; (2) the inversion algorithm to iteratively calculate the voxel intensity based on PSF and projections from all perspectives. In FCICT, the inversion process uses the algebraic reconstruction technique (ART) which has been extensively validated by our past works (Ling et al. 2021(Ling et al. , 2020. Hence, the following describes the establishment of PSF matrix for confined-space tomography. The demonstration of PSF calculation is shown in Fig. 1. There are 2 steps of the calculation process: first, the reverse ray-tracing to establish the reverse mapping relationship from 2D pixel center points to the 3D measurement domain following the Snell's Law; second, linear interpolation in the 3D domain to obtain the forward mapping relationship (i.e., the ray-tracing) from 3D voxel centers to corresponding projections. It's worth to be noticing that since the parameters of cylinder have been considered in the process of ray-tracing, the thickness of cylinder or other parameters have no effect on the calculation accuracy of the optical path. Figure 1a first demonstrates the process of reverse ray-tracing. The region surrounded by red dotted line is a partial enlarged schematic of a set of adjacent 4-pixel which's center points are I′, J′, K′ and L′, respectively. Consequently, the rays R I , R J , R K and R L radiated from I′, J′, K′ and L′ pass through the center of lens (pinhole model) and then are refracted consecutively by the outer and inner walls of optical cylinder accorded with Snell's law. The detailed calculation of refracted ray directions based on Snell's law can be found in our previous work (Ling et al. 2020). Point M is an arbitrary point in measurement domain surrounded by the adjacent 4-rays. The distances from M to 4-rays are calculated as MI, MJ, MK and ML, respectively. Hence, the distances between M and IJ, JK, KL and LI are l 1 , l 2 , l 3 and l 4 , respectively. As the next step, the projecting point of M on camera chip (i.e., M′) can be calculated by the following linear interpolation: where l 1 ′, l 2 ′, l 3 ′ and l 4 ′ are the distance from point M′ to I′J′, J′K′, K′L′ and L′I′, respectively. Successively, the ray-tracing relationship between point M and its projecting point M′ is established. Following similar process, the reversed ray-tracings are calculated from all adjacent 4-pixel groups on the camera chip, then the one-to-one correspondence from any point in measurement domain to its projection on camera chip is mathematically determined. The next step is to decide the projection region of voxels according to above one-to-one correspondence. As shown in Fig. 1b, points ABCDEFGH are the corner points of a voxel centered at point N. Their projections on camera chip are calculated to be points A′B′C′D′ E′F′G′H′, respectively. Thus, the projection region of voxel N is decided by the maximum area covered by above projection points (e.g., E′F′G′H′ in Fig. 1b). The intensity P of a pixel contributed by voxel N is then calculated by: where V is the voxel intensity, D lens and f are the diameter and f-number of the lens, respectively. l represents the distance from voxel center point N to the lens center. S cover is the Demonstration of the FCICT imaging algorithm. a Reversed ray-tracing process: The rays emitted from adjacent 4-pixel center points are respectively R I , R J , R K and R L . M is an arbitrary point in the measurement domain surrounded by the four rays, and its projection point is M′ calculated by linear interpolation; b The ray-tracing from voxel ABCDEFGH to its projection A′B′C′D′E′F′G′H′ which presents the signal mapping relationship from a voxel to pixels area of pixel covered by the projection of voxel N, and S pixel is the total area of a pixel. Following Eq.
(2), the mapping relationships from voxels to pixels are established, and the voxel intensity V is then iteratively solved by ART. The effects of cylinder thickness on the reconstruction quality have been quantified accurately in the FCICT algorithm, and can be neglected in practical cases (under the assumption that the refraction index is uniform). The effectiveness of FCICT imaging algorithm and the accuracy of the experimental setup are validated by a proof-of-concept calibration plate imaging experiment, as shown in Fig. 2. The calibration plate is installed at the meridian plane of an optical cylinder. The plate has chessboard pattern and consists of black-and-white rectangular grids with a number distribution of 28(horizontal) × 17(vertical). The size of each grid is 3.0 mm × 3.0 mm, leading to a total area of the plate of 84 mm × 51 mm. The cylinder is made of quartz with uniform refractive index of 1.4, and has a size of 131 mm (outer diameter) × 92 mm (inner diameter) × 42 mm (height). In Fig. 2a, a projection of a calibration plate is captured, and a portion of it is covered by the cylinder while the rest is in open-space. The focus point in Fig. 2a is the red point h at the center of image. The projection covered by red line is the distortion area affected by the cylinder. The blue vertical auxiliary lines mark the leftmost distorted projection, center line and rightmost undistorted projection, respectively. The orange horizontal auxiliary lines mark the center horizontal line and a quasi-horizontal line at the top of the last row of grids in the projection, respectively. The vertical distances between their left and right ends are d 1 and d 2 , respectively. Figure 2b is the simulated projection at the same observation orientation as Fig. 2a. According to visual comparison of vertical auxiliary lines, the measured and simulated projections match well both in undistorted and distorted area. The orange auxiliary lines also mark the same location of calibration plate, and the corresponding distance are d 1 ′ and d 2 ′, respectively. Compared d 1 ′ to d 1 and d 2 ′ to d 2 , the averaged absolute error is less than 0.2 pixels, and the relative error is less than 0.5%, normalized by the size of grid (i.e., 40 pixels). Through the above experiments, the relationship between distortion and observation orientation is also obtained. The horizontal distortion would be large on the projection region close to the cylinder inner wall, while the vertical distortion would significant increase if the detector is located with large inclination angles. More detailed analysis can be found in the ref (Ling et al. 2020). By now, the effectiveness of FCICT imaging algorithm and the accuracy of the experimental setup are validated using a static phantom.
The FCICT method is numerically validated using a simulated 3D phantom. The measurement domain of the numerical validation is decided in a Cartesian coordinate system with a dimension of 92 mm (X direction) × 92 mm (Y direction) × 40 mm (Z direction), sufficient to encompass the area inside the optical cylinder. The measurement domain is discretized into 184 voxel × 184 voxel × 80 voxel with the nominal resolution of 0.5 mm/ voxel. The 3D phantom consists of 2 axisymmetric hollow objects as shown in Fig. 3a. The height of both objects is 60 mm. The vertical axes are located at X = 15 mm, Y = − 20 mm and X = − 15 mm and Y = 20 mm, respectively. The phantom intensity on the surface is set as 10 (arbitrary value), and all internal voxels are set as 0. Figure 3b-i present a group of 8 simulated projections of the phantom by the FCICT imaging model considering the distortion of the cylinder refraction. The magnification of all views is set to 11. The azimuth angles of 8 views are evenly distributed in the range from 0° to 180° with a uniform increment of 25.7°, and the inclination angles are set to 0° for all views. A random angle error within ± 0.1° is added to the azimuth angles to simulate the experimental uncertainty. All projections have a pixel resolution of 420 pixel × 180 pixel which is sufficient to enclose the entire measurement domain. Figure 4a shows the FCICT tomographic reconstruction of the phantom using the projections in Fig. 3b-i. As presented in Fig. 4a, the overall structure of phantom is promisingly reconstructed, though the objects present slightly rough surface (as also indicated by the projections in Fig. 3b-i caused by the manually added random error. To better express the internal details, Fig. 4b, c show the 2D phantom slices cropped at Y = − 20 mm and Y = 20 mm, respectively. As can be observed, the phantom boundary signal is shape and clean (with acceptable uniformity). As a comparison, Fig. 4d shows the reconstructed phantom processed by traditional open space CT algorithm (still using distorted projections as input). It can be observed that there is significant noise at the surface front of both objects. Furthermore, Fig. 4e, f show the 2D slices at Y = − 20 mm and Y = 20 mm, respectively. It is also clearly observed that the surface front has noticeable inhomogeneous intensity, contradicted to the given uniform phantom distribution. For quantitatively comparing the quality of reconstruction results, a correlation coefficient developed from (Krijnen 1994) is used to express the similarity between reconstructed results and original 3D phantom as shown in Eq. (3): where V phan and V recon are the 3D signal matrix of the phantom and reconstructed flames, respectively; the overline above a matrix indicates the average of all non-zero values of the matrix. The closer Corr is to 100%, the higher similarity is between the phantom and the reconstruction. The correlation coefficient is 88.30% for the FCICT reconstruction and

Experiment Setup and FCICT Reconstruction
This section reports the experiment setup and tomographic reconstruction using FCICT on a two-phase flow in confined-space (i.e., inside an optical engine cylinder). The schematic setup of the experiment is shown in Fig. 5. All optical instruments are installed on an optical platform. The optical cylinder is set at the center of optical platform held by a support seat. A water liquid jet (atomized by ultrasonic device, model: NRWT-1) is set inside the optical cylinder. The frequency and power of the ultrasonic device is 1.7 MHz and 4 W, respectively. The physical dimension of atomized droplets is about 5 μm. The height of liquid jet output is set the same as the bottom of optical cylinder. The horizontal location of liquid phase can be adjusted free in the optical cylinder. The liquid jet can be approximately regarded as a reflecting surface, since the aim of this research is the application of new method in measuring the 3D surface properties of the two-phase flow. Hence, both multi-scattering and absorption among liquid droplets are neglected, since they are out of the scope of the paper. The assessment of developed CT method in more complex flow fields is scheduled in our future works. A continuous wave laser (RayPower Fig. 5 Experiment setup. A liquid jet is placed inside an optical cylinder mounted on an optical platform, surrounded by a circular rail. In total of 8 fiber detectors on 2 four-to-one fiber bundles are mounted on the circular rail to transmit the liquid jet projections illuminated by a continuous-wave laser to the fiber outputs. 2 synchronized high speed cameras are used to capture the liquid jet projections from the fiber outputs, and store the projections in a lab computer 5000, wavelength 532 nm) is set on the optical platform to illuminate the liquid jet signal. The ray output from laser is reflected 3 times by a series of 45° reflecting mirrors, then expanded consecutively by a couple of plano-concave lenses (f = 50 mm and 150 mm) into a conical beam to illuminate the entire measurement field inside the cylinder. A circular rail with the diameter of 1200 mm is set coaxial with the optical cylinder. There are in total of 8 fiber detectors installed on the circular rail, all equipped with Nikon MF 50 mm/f 1.4 lens, to record projections from different orientations simultaneously. The fiber detectors are connected to a couple of four-to-one fiber bundles, so the output images ( Based on above experimental setup, the FCICT algorithm is utilized to reconstruct the 3D flow patterns of the two-phase flow at different time sequences. The camera repetition rate is set to 300 Hz with the exposure time of 3.33 ms. The power of laser is set to 1.5 W. During the experiments, the process of the two-phase flow emerging from the jet outlet is performed and recorded repeatedly, resulting in a total of 21 selected sets of projections, each with 60 frames to enclose the flow development. Figure 6 shows one frame of a set of projections for presentation. The recorded liquid jet projections have uniform pixel dimensions of 270 pixel × 220 pixel, cropped from the original camera chip. The azimuth angles (i.e., the angle between the projection of optical axis of camera on O-XY plane and X axis) of projections recorded through Fiber 1 to 8 are 9.6°, 34.3°, 60.9°, 84.2°, 105.7°, 132.7°, 157.7° and 183.6°, respectively. The inclination angles (i.e., the angle between optical axis of camera and O-XZ plane) are − 2.9°, − 2.5°, − 1.0°, − 1.2°, 0.4°, 3.8°, − 2.7° and − 0.7°, respectively. The azimuth and  (Kang et al. 2014) that has been extensively applied in our past works. As shown in Fig. 6, the yellow line represents the location of the optical cylinder top surface where the signal is ignored. This is because the optical cylinder has chamfers at the junction of the top and bottom surfaces to avoid stress concentration problems. The chamfer is a frosted surface, so the direction of light propagation though it can not be predicted. This region (with ~ 15 pixels in height) is expelled from reconstruction and further discussion. Comparing the projection signal distribution inside and outside the cylinder, it is noticeable that the liquid jet signal intensity outside the cylinder is evidently stronger than that inside the cylinder, since the ray-tracing is partially absorbed or reflected by the cylinder. Moreover, the liquid jet generally has straight-up pattern inside the cylinder, but tends to deform in the horizontal directions when coming out of the cylinder. One possible reason is that the outer part of the liquid jet flow has enhanced diffusivity in the open space with higher flow unsteadiness, though we do not intendedly generate disturbance in the ambience.
Based on a set of liquid jet projections in Fig. 6, the 3D reconstruction is performed, as shown in Fig. 7. The 3D measurement domain is 92 mm × 92 mm × 78 mm and discretized into 184 voxel × 184 voxel × 156 voxel, with the nominal spatial resolution 0.5 mm. The height of measurement domain is set larger than the simulation case in Sect. 3 to enclose the portions of liquid jet outside the cylinder. As shown in the 3D renderings in Fig. 7, the regions circled out by yellow present the sub-domain inside the cylinder. In Fig. 7a, the sub-domain inside the cylinder is reconstructed by FCICT and it outside cylinder is reconstructed by traditional open-space computed tomography (named CT later for brevity). The red curve in Fig. 7a marked out parts of the liquid jet outline (viewed from the azimuth angle of 85°). It is corroborated by the curve that the inner portion of the liquid jet spray has smooth transition to its outer portion. As a contrast, in Fig. 7b, the measurement sub-domains both inside and outside the cylinder are reconstructed by CT. There is significant signal discontinuity and pattern displacement between the 3D renderings inside and outside the cylinder, which is contradicted to the flow continuity. Besides, increasing artificial noise can be found in the 3D reconstruction inside the cylinder, since CT fails to correct the refracted ray-tracings from voxels to pixels, leading to intensity mismatch from different projections and eventually Fig. 7 3D liquid jet structure reconstructed by a FCICT, with continuous and smooth intensity distribution; b CT with discontinuous intensity between sub-domains inside and outside the cylinder to increased reconstruction error. In short, the comparison between Fig. 7a, b provides qualitative validation for the correctness and necessity of FCICT.

Surface Curvature and Velocity Measurements
On the basis of sequentially reconstructed 3D liquid jet structure, the 3D curvatures and velocity distributions on the surface of the two-phase flow are calculated and demonstrated in this section. The curvature of liquid jet surface determines the surface tension force according to the Laplace equation. Besides, the surface tension force is determined by curvature gradient. In atomization systems, the liquid surface is the balance of inertia and surface tension, as mathematically represented by Weber number. From the past researches, the surface tension has a significant influence on the fluid flow with free surface or multiphase flow with sharp interface (Zhang 2010). In past 2D flow or flame studies, 2D curvatures were first extracted on captured images, then proved to be correlated with the flow or flame essence. For example, the flame curvature influences the local heat release rate (Kosaka et al. 2020) and flame displacement speed (Sinibaldi et al. 1998;Tsuchimoto et al. 2009). However, since most practical flows or flames are inherently 3D with asymmetric surface distribution, the 3D curvature measurement is a necessity to characterize the flame nature. Recently, multiple calculation methods of 3D curvature were investigated and developed, such as the central difference method in quad-plane PIV measurements (Kerl et al. 2013), the three-point finite-differencing scheme in (Ma et al. 2016b;Yu et al. 2020), the local polynomial fitting method in flame surface measurements (Wiseman et al. 2017), and so on. By implementing these methods, 3D curvatures were calculated and discussed in some studies to reveal interesting correlations with the flow or flame essence. For example, Chi et al. finds that the response of 3D flame propagation speed to the 3D mean curvature is positively correlated (Chi et al. 2022). Enlighted by past efforts, this work investigates the 3D curvature distribution on the surface of reconstructed two-phase jet flow, aiming at finding quantitative relationship with the flow dynamics, represented by the surface velocity so as to show the potential of the FCICT in 3D measurement and 3D analysis of practical optical signals.

3D Curvatures on Two-Phase Flow Surface
Two different 3D curvatures, mean curvature (K Mean ) and Gaussian curvature (K Gauss ), are calculated based on the tomographic reconstructions of the liquid jet surface found in Sect. 3. For arbitrary point on the liquid jet surface, infinite number of 2D curvatures can be decided by slicing the surface with a 2D plane at the location of the point. Among all 2D curvatures, the maximum and the minimum are named principal curvatures. The average and the product of two principal curvatures are called mean curvature (K Mean ) and Gaussian curvature (K Gauss ), respectively. Both K Mean and K Gauss can be only obtained by 3D measurements, and are usually used to describe the instantaneous surface shape and forecast the flow development (Chi et al. 2022). The process of K Mean and K Gauss calculation is divided into 4 steps. First, a point cloud of liquid jet surface is extracted according to the tomographic reconstruction. Second, the normal vectors of all points included in the point cloud are respectively calculated. Third, the first and the second fundamental forms of the surface fitted by the point cloud are calculated 1 3 through the normal vectors obtained in the second step. Finally, the distributions of K Mean and K Gauss are derived from the first and the second fundamental forms.
Specifically, in the first step the point cloud is extracted from the 3D reconstruction using the ISO-surface function in MATLAB, as shown in Fig. 8. The point cloud denotes all surface points of the reconstructed liquid jet. For an arbitrary point P in the cloud, the closest n points around P (namely Q 1 , Q 2 , …, Q n ) are marked out, as shown in Fig. 9. Hence, the curved surface s is fitted through points Q 1 -Q n to enclose point P, expressed as Eq. 4: where u and v are coordinates on the point cloud, which are a set of linearly independent vectors. Equation 4 determines the intrinsic geometry of the surface at P. Moreover, the derivative of s is shown in Eq. 5: where r u and r v are the differential coefficient of s on coordinate u and v, respectively. Hence, the coefficients in the first fundamental forms E, F and G of the surface s can be calculated by Eq. 6: Besides, the normal vector n of surface s is obtained by Eq. 7: The second fundamental forms L, M and N can be expressed by Eq. 8: where r uu represents the partial derivative of r u on u coordinate, similarly for r uv and r vv . According to Eq. 6 and 8, the K Mean and K Gauss can be calculated by Eq. 9-10, respectively (Chern 1945): Mean and Gaussian curvatures are used to distinguish local features of the surface. By calculating the mean curvature, a surface point can be classified as a concave point (K Mean < 0) or a convex point (K Mean > 0). Similarly for Gaussian curvature, a surface point can be defined as elliptic point (K Gauss > 0) or a hyperbolic point (K Gauss < 0). Combining both curvatures, parabolic points (K Gauss = 0 and K Mean ≠ 0) and planar points (K Gauss = 0 and K Mean = 0) are defined. Figure 10 presents the K Mean and K Gauss distributions of the reconstructed twophase flow in Fig. 8. As shown in Fig. 10a, the portions of liquid jet outside the cylinder (47 mm < Z < 60 mm, in red) and inside the cylinder (15 mm < Z < 45 mm, in purple) are respectively studied. The curvature results between the studied regions are corrupted by the cylinder edge, thus expelled from the discussion. Generally, the mean curvatures on the outside liquid jet portion have moderate value (− 0.14 < K Mean < 0.16) compared to those on the inside portion (− 0.29 < K Mean < 0.3). One possible reason is that kinetic energy of liquid jet keeps reducing it emerges from the jet outlet. After the vaper develops outside the cylinder, it has experienced sufficient expansion and tends to homogenously distributed in space, leading to even reduced absolute value of K Mean (i.e., the surface tends to be flat). On the contrary, higher flow velocity and more robust air entrainment tend to generate wrinkled surface on the two-phase flow inside the cylinder, leading to noticeably increased K Mean . In this comparison, it is noteworthy that the sign of K Mean is related with the definition of surface coordinates, so the absolute value of K Mean is more concerned. Similar observation can be found in Fig. 10b for the Gaussian curvatures. A quasi-uniform K Gauss distribution is observed on the outer liquid jet portion, with very slight variation of − 0.005 < K Gauss < 0.0015. Different from mean curvature, Gaussian curvature demonstrates the intrinsic character of surface. Therefore, the sign of K Gauss is independent from the surface coordinates. That is to say, the points on the outer surface are primarily elliptic.
Such distribution is in accordance with the Fick's law validated by past works (Fick 1855;Porteous 2001). Contrarily, the Gaussian curvatures inside the cylinder are distributed in wider range (− 0.05 < K Gauss < 0.047) due to strong mixing effect between the air and water droplets. Such impact also causes elliptic and parabolic points alternately distributed in a narrow region (70 mm < X < 80 mm, 30 mm < Y < 35 mm, 23 mm < Z < 33 mm, as marked out in Fig. 10b). The statistics of 3D curvature distributions are analyzed based on all 21 sets of liquid jet frames recorded at 300 Hz (60 frames per set). Figure 11 presents the curvature statistical Fig. 10 The a mean and b Gaussian curvature distributions on the liquid jet surface Fig. 11 The variations of a the numbers and b the probability distribution of surface points over time results for the set that the liquid jet in Fig. 8 belongs to. As shown in Fig. 11, the total recording period for liquid jet development is 200 ms, where the origin (0 ms) is set as the frame in which the liquid jet just leaves the jet. Figure 11a shows the population of classified surface points. As for the concave and convex points (distinguished by K Mean ), both numbers significantly increase with time during the liquid jet development. This is because the liquid jet surface area is expanded during its upward movement. It is also observed that the populations of concave and convex points grow faster after ~ 80 ms, when the liquid jet head leaves the cylinder. The higher growth rate is supposed to be dominated by enhanced liquid jet diffusion outside the cylinder. Similar trends are also observed on the population variations of the elliptic and hyperbolic points (distinguished by K Gauss ). More elliptic points than hyperbolic points are observed at all time sequences, while the number difference tends to decrease over time. This may be caused by the liquid jet diffusion process with irregular wrinkles generation on the surface (the readers are referred to Video 1 for straightforward observation). Specially, no parabolic or planar point (distinguished by K Mean and K Gauss ) is obtained among all tested cases. To more distinctly compare the population variation of different surface points, Fig. 11b shows the development of probability distribution of all points, calculated based on Fig. 10a. Specifically, both proportions of concave and convex points fluctuate around 50%, since the liquid jet is in the state of Brownian Motion (Feynman 1964;Chern 1945). Comparatively, the elliptic points occupy ~ 90% of all surface points at 0 ms whereas the hyperbolic points only occupy 10%. This is because the liquid jet head diffuses upward in the shape of quasi-hemisphere when it just leaves the jet, leading to majority of elliptic points. From 0 to 200 ms, the proportion of elliptic points gradually decreases from 90 to 50%, since the initial kinetic energy is exhausted over time while the irregular Brownian Motion becomes dominant factor in the liquid jet diffusion process (Feynman 1964).

Liquid jet Surface Velocity
The liquid jet surface velocity distribution is calculated using so-called normal-vector method (Wiseman et al. 2017) as illustrated in the 2D schematic plot Fig. 12. As shown in Fig. 12, surface s 1 and s 2 are consecutive frames of a liquid jet surface (divided by unity time). For s 1 , the direction of surface normal vector is defined from inside (the shaded area surrounded by red line) to outside of the liquid jet. Specifically, for point P 1 on surface s 1 , the normal vector is denoted by n 1 , which then intersects surface s 2 at P 1 ′. Hence, the surface velocity vector on point P 1 can be defined by V 1 (green vector). The velocity is positive since the n 1 and V 1 are in the identical direction. In contrast, point P 2 has its normal vector n 2 and velocity vector V 2 in opposite directions, leading to negative velocity value. Although the normal-vector method is limited in measuring Fig. 12 The normal-vector method for calculating the surface velocity distribution velocity of highly wrinkled surface (Chi et al. 2022), the normal-vector method in this work is considered to be reliable, because the jet flow in this work is measured under quiescent ambience, resulting in the absolute value of the mean curvature less than 0.3 for all cases.
By implementing above method, the surface velocity distributions of the liquid jet are calculated. Figure 13a shows the surface velocity distribution of the liquid jet reconstructed in Figs. 6,7 and 8 (recorded at 116.55 ms). In details, for the portion outside cylinder in Fig. 13a, the velocity of the liquid jet circled by green triangle is mainly negative except the region in the lower right corner. Meanwhile, the velocities in other regions are mostly positive. For the portion inside cylinder, the velocity distribution is mostly negative or close to zero, as represented by the green rectangle. To better demonstrate the deformation trend of liquid jet surface as a consequence of surface velocity, Fig. 13b superimposes the reconstructed liquid jet surface in adjacent frames (i.e., 116.55 ms in blue and 119.88 ms in red). By comparing Fig. 13a and b, the surface velocity distribution overall matches the surface location reconstruction (e.g., similarity distribution marked out by green triangle and rectangle). Specifically, 99.7% of surface points with positive velocity in Fig. 13a are also painted in red in Fig. 13b, while 99.5% of negative velocity points are presented in blue. Such comparison indicates the consistency between the normal-vector method and the FCICT method and good robustness of normal-vector method.
After the 3D curvatures and surface velocity distributions are determined, their relationship is investigated. Figure 14 shows the correspondence between the average surface velocity (absolute value) and mean/Gaussian curvatures by examining all surface points of reconstructed liquid jet at 116.55 ms. Specifically, in Fig. 14a, the horizontal axis (K Mean ) is ranged from − 0.25 to 0.25 with an increment of 0.01. The vertical axis presents the average velocity of surface points whose mean curvature is contained in each K Mean element. Such statistical processing is performed to reduce the influences by extreme curvature/velocity values and to reveal the major relationship. It is noticeable where |V| represents the average absolute surface velocity. Similar statistical processing is also performed for Gaussian curvature-surface velocity relationship shown in Fig. 14b, except that the horizontal axis (K Gauss ) is ranged from − 0.003 to 0.003. As illustrated by Fig. 14b, |V| presents monotonical increasing trend with larger K Gauss for elliptic points (K Gauss > 0). Corresponding numerical fitting is presented in Eq. (12). Besides, it is also noticed that no obvious correspondence is observed for hyperbolic points.
To further understand the correspondence of surface velocity and 3D curvatures in Fig. 14, we propose a qualitative explanation by considering the movement directions of different surface points of the liquid jet. In brief, a surface point with higher level of concavity or convexity would be more likely to possess higher velocity. As shown in the schematic plot in Fig. 15, point S 1 is assumed to obtain high concavity or convexity during time t 0 → t 1 . At t 0 , S 1 coincides with location S 0 , and the region around S 0 has very small mean curvatures (simplify to 2D, represented as the green line). n is the normal vector from point S 0 . Two situations are likely to occur to form S 1 : (1) the displacement (S 0 → S 1 ) is neglected whereas all surrounded points simultaneously move in the opposite direction of n (Fig. 15a); (2) the displacement (S 0 → S 1 ) is large while surrounded points move slightly (Fig. 15b). Since the surface points, initially driven by the kinetic energy from the jet outlet, are increasingly affected by the Brownian Motion when the liquid jet gradually develops into the open space, the possibility would be reduced for a group of close points to move in the same direction. On the contrary, the possibility would increase for a single point to gain excessive velocity and leave surrounding points behind. If that is the case, the second situation (Fig. 15b) would be more likely to occur. In 3D space, increased 2D curvatures of point S 1 would be consequently obtained from all directions, leading to larger K Mean . Moreover, all 2D curvatures on S 1 are likely to have the same sign, indicating that S 1 is an elliptic point. This assumption is in accordance with the observation in Fig. 14b that only elliptic points are correlated with K Gauss . Furthermore, statistical analysis is performed to quantify the correspondence between surface velocity and 3D curvatures. To aid the analysis, we denote the maximum absolute surface velocity, absolute mean curvature and Gaussian curvature as |V max |, |K M,max | and K G , max . By analyzing all surface points on studied cases, it is found that about five-sixths of surface points with |K Mean |> 0.9|K M,max | or K Gauss > 0.8 KG, max have the surface velocity |V|> 0.7|V max |. Such result further proves the correspondence between V, K Mean and K Gauss , all calculated based on FCICT tomographic reconstructions.

Summary and Conclusion
In summary, this work reports a full-field cross-interface tomography algorithm (FCICT), and the emphasis on its numerical validation and practical applications. The FCICT utilizes the Snell's law and reverse ray-tracing to obtain the mapping relationship between 2D projections and 3D optical field under the impact of imaging distortion caused by an optical engine cylinder. Sequentially, the ART method is used to iteratively solve the 3D optical field after view registration. In the numerical validation, a 3D phantom is established with its projections evenly distributed from 0° to 180° (azimuth angle). The simulated projections are then input into the FCICT reconstruction algorithm with ± 0.1° angle error added. The FCICT reconstruction shows similar signal distribution as the original phantom, with high correlation coefficient of 88.30%. In contrast, the reconstruction by traditional CT shows coarse reconstruction, with artificial noise, pattern displacement and low correlation coefficient of 37.06%. After numerical validation, the FCICT is then employed in practical reconstructions of a two-phase flow generated from the optical cylinder. A water liquid jet is generated inside the cylinder and gradually move upwards to leave the cylinder. Both liquid jet portions inside and outside the cylinder are illuminated by a CW laser. The Mie scattering signal is simultaneously captured by 8 fiber detectors and recorded by 2 high speed cameras through 2 four-to-one fiber bundles. By comparing the flow continuity and signal intensity distribution in the reconstructions, the FCICT shows better performance The results show that the population of concave/convex points (denoted by mean curvature) and elliptic/hyperbolic points (by Gaussian curvature) all increase with time, since the liquid jet surface area is expanded during its upward movement. Besides, with irregular wrinkles appear on the liquid jet surface as it develops outside the cylinder, the population of hyperbolic surface points increases while that of the elliptic points reduces. The probability distribution analysis shows that the proportions of concave and convex points are both about 50% as the liquid jet is fully developed, similar to those of elliptic and hyperbolic points. Such phenomenon can possibly be attributed to the gradual dominance of Brownian motion. Finally, the relationship between surface velocity and 3D curvatures is investigated, showing that the absolute surface velocity monotonically increases with the absolute K Mean and K Gauss (for elliptic points only). Such phenomenon can also be attributed to the dominant Brownian motion that empowers random surface points to move away from nearby points and turn into concave/convex points in all 2D directions (i.e., elliptic points). In sum, the FCICT algorithm developed in this work is proved to be capable of solving confined-space tomography problems. During 3D measurements of a two-phase jet flow, the 3D curvature and velocity distributions are achieved from the accurate topology reconstruction by FCICT, providing quantitative evaluation of the flow field.