An improved MLEM method applied to the reconstruction of the three-dimensional bubble flow field

Abstract

The present invention belongs to the field of three-dimensional flow field reconstruction of bubbles, and relates to an improved MLEM method applied to the three-dimensional flow field reconstruction of bubbles. This method combines multi-view vision measurement with the improved MLEM method. First, multiple sets of calibration plate images captured synchronously by multiple cameras are used to calibrate the multi-camera system, and the world coordinate systems of the multiple cameras are unified. Then, multiple cameras are used to synchronously capture the bubble flow field to obtain bubble images at multiple angles at the same moment. Finally, in view of the characteristics of backlight illumination of the bubble flow field, the method for setting the initial iteration value and the projection and back-projection methods in the iterative process in the MLEM method are improved, realizing the rapid reconstruction of the three-dimensional flow field of bubbles based on the improved MLEM method. The improved MLEM method designed by the present invention for the three-dimensional flow field reconstruction of bubbles not only improves the calculation speed but also enhances the reconstruction accuracy, and can be widely applied to the field of three-dimensional flow field reconstruction of bubbles in gas-liquid two-phase flows.

Description

Technical Field

[0001] The present invention relates to an improved MLEM method applied to the three-dimensional flow field reconstruction of bubbles. More specifically, the present invention relates to a method for reconstructing the three-dimensional flow field of bubbles based on the combination of multi-view vision measurement and the improved MLEM method. Background Art

[0002] Gas-liquid two-phase flow, as one of the most common and prevalent forms of two-phase flow, widely exists in human production, life, and practical activities. Bubbles, as an important part of the gas-liquid two-phase flow field, their characteristics such as shape and distribution are closely related to the efficiency of industrial production. Studying the reconstruction technology of the three-dimensional flow field of bubbles is of great significance for studying the motion characteristics of bubbles, mastering the flow mechanism of gas-liquid two-phase flow, and promoting the development of industrial production.

[0003] The three-dimensional reconstruction technology based on images has been widely used in fields such as medical measurement and space flight due to its advantages such as intuitive measurement and low cost. The iterative reconstruction algorithm is a reconstruction method applied to computer tomography. Using the tomography method, a three-dimensional reconstruction problem can be simplified into a two-dimensional reconstruction dimensionality reduction model. In the field of flow field reconstruction, the iterative reconstruction algorithm has also been applied to a certain extent. The traditional iterative reconstruction algorithm is mainly based on the projection value of light intensity superposition, and the light intensity distribution of voxels in space is inversely calculated through iteration. However, the traditional iterative reconstruction algorithm has strict requirements for the projection angle during the three-dimensional reconstruction of images, the equipment cost is high, and the reconstruction accuracy and speed are not ideal. Moreover, during the image acquisition process of gas-liquid two-phase bubbly flow, the backlight illumination method is usually used to capture the characteristics of bubbles. When the ray passes through the flow field and encounters the target bubble, a corresponding projection point will be formed on the camera, rather than the direct superposition of light intensity. At this time, the traditional iterative reconstruction algorithm is not fully applicable and needs to be improved to achieve the three-dimensional reconstruction of the gas phase distribution in the bubble flow field. Summary of the Invention

[0004] In view of the above deficiencies, the present invention designs an improved Maximum Likelihood Expectation Maximization (MLEM) method applied to the three-dimensional flow field reconstruction of bubbles. First, calibrate the multi-camera system to unify the world coordinate systems of the multi-cameras; then, aiming at the characteristics of backlight illumination of the three-dimensional flow field of bubbles, improve the setting method of the iterative initial value and the projection and back-projection methods in the iterative process of the MLEM method. Finally, an improved MLEM method applied to the three-dimensional flow field reconstruction of bubbles is designed, realizing the rapid reconstruction of the phase distribution of the three-dimensional flow field of bubbles.

[0005] The hardware system of the improved MLEM method applied to the three-dimensional flow field reconstruction of bubbles includes: K high-speed CCD cameras for image acquisition, with a resolution of 1024×1024 pixel and a shutter speed of 1 / 2000 s, and the K cameras shoot synchronously; 1 computer for saving data images; 1 checkerboard calibration board for camera calibration; 1 air pump for generating bubbles; and 1 water tank for simulating gas-liquid two-phase flow.

[0006] The present invention designs an improved MLEM method applied to the three-dimensional flow field reconstruction of bubbles, which is characterized by including the following steps:

[0007] Step 1: Place the K cameras around the shooting target and on the same horizontal plane, and the overall shooting angle of the cameras is less than 180 degrees.

[0008] Step 2: Place the checkerboard calibration board in the target area to be photographed, change the position of the checkerboard calibration board multiple times, and use the K cameras to synchronously photograph the checkerboard calibration board to obtain multiple groups of checkerboard calibration board images at K angles at the same moment.

[0009] Step 3: Perform single-object calibration on each camera to obtain the internal and external parameter information of each camera.

[0010] During calibration, the computer will obtain the corner point information of the checkerboard images taken by each camera, and each camera defines the coordinates of the checkerboard corner point closest to the upper left corner as the origin of the world coordinate system, the horizontal direction of the checkerboard as the X-axis, the vertical direction of the checkerboard as the Y-axis, and the direction perpendicular to the checkerboard plane as the Z-axis. Taking this as the world coordinate system, the external parameters will be calculated based on this world coordinate system. Since the position of the checkerboard needs to be changed during calibration and multiple images are used for calibration, the world coordinates corresponding to different images are also different, and the calculated external parameters are also different. In order to make the world coordinate systems of the K cameras have a unified reference, the external parameters of each camera are calculated based on the calibration images synchronously taken by the K cameras at a certain moment. The internal parameters of the camera represent the relationship between the camera coordinates and the pixel coordinates and are only determined by the camera itself and can be uniquely determined.

[0011] Perform single-object calibration on each camera respectively to obtain the internal parameter matrix M n 、rotation matrix R n and translation vector T n of the nth (1 <= n <= K) camera. At this time, the rotation matrix and translation vector of each camera are both relative to the world coordinate system determined by the checkerboard calibration board at the same position.

[0012] Step 4: Use camera 1 as the reference camera, and use the internal and external parameters of the cameras obtained in step 3 to calculate the position relationship parameters between other cameras and camera 1.

[0013] Perform binocular calibration on Camera 1 and Camera n, and calculate the position relationship parameters of Camera n relative to Camera 1; assume that a certain point in space is q, and its world coordinates are known as q w , and the coordinates of point q in the camera coordinate systems of Camera 1 and Camera n are

[0014] q 1 = R 1 ·q w + T 1 Equation (1)

[0015] q n = R n ·q w + T n Equation (2)

[0016] From Equation (1) and Equation (2), the relationship between the two camera coordinates q 1 and q n is

[0017] q n = R 1n ·q 1 + T 1n Equation (3)

[0018] where

[0019]

[0020] T 1n = T n - R 1n ·T 1 Equation (5)

[0021] R 1n and T 1n are respectively the rotation matrix and translation vector between Camera 1 and Camera n, which are the position relationship parameters between cameras to be found.

[0022] Step 5: Based on the world coordinate system of Camera 1, use the position relationship matrix between Camera 1 and Camera n described in Step 4 to convert the world coordinate system of Camera n into the world coordinate system of Camera 1, then the external parameter matrix of Camera n is

[0023]

[0024] Step 6: The relationship between the world coordinate system (X w , Y w , Z w ) and the pixel coordinate system (u, v) is shown in Equation (7).

[0025]

[0026] Among them, Z c is the coordinate value of the current coordinate point on the Z-axis of the camera coordinate system; M is the internal parameter matrix of the camera, including f x , f y , u 0 , v 0 , s and other parameters, f x and f y represent the scale factors for the transformation between the pixel coordinate system and the camera coordinate system of the camera, (u 0 , v 0 ) is the corresponding coordinate of the origin of the image coordinate system in the pixel coordinate system, and s is the skew factor; N is the external parameter matrix of the camera, including the rotation matrix R and the translation vector T; MN is the camera projection matrix, and the matrix dimension is 3×4.

[0027] Substituting the internal parameters of each camera calculated in step 3 and the external parameter matrix of each camera calculated in step 4 into formula (7), the projection relationship expression between the world coordinate system and the pixel coordinate system of each camera can be obtained, as shown in formula (8).

[0028]

[0029] In the formula, (u n , v n ) is the pixel point coordinate in the pixel coordinate system of camera n; M n is the internal parameter matrix of camera n.

[0030] Step 7: After completing the content described in step 6, keep the relative positions of the cameras unchanged, place a water tank in the photographed area, introduce bubbles into the water tank, and use K cameras to synchronously photograph the bubble flow field to obtain bubble images at K angles at the same moment.

[0031] Step 8: Based on the K-angle bubble images collected in step 7, perform denoising and binarization processing on the images to obtain the binarized bubble images I 1 , I 2 ,..., I n ,..., I K .

[0032] Step 9: Set the maximum number of iterations Ite m and the minimum allowable error error min .

[0033] Step 10: Represent the projection matrix MN in formula (7) described in step 6 with formula (9).

[0034]

[0035] Among them, C ab is the matrix element in the projection matrix, a = 1, 2, 3, b = 1, 2, 3, 4.

[0036] Then, the projection relationship from the world coordinate system of the camera to the pixel coordinate system can be expressed in the form shown in formula (10).

[0037]

[0038] Stratify the three-dimensional space region along the Z w axis of the world coordinate system. Consider the three-dimensional space region as a continuous multi-layer two-dimensional image. For the known Z w layer, its back-projection calculation formula is as shown in formula (11).

[0039]

[0040] Substitute the two-dimensional pixel coordinates (u, v) and the Z w coordinate of the current layer into the back-projection calculation formula, and the X w and Y w coordinates of the corresponding point on the Z w layer of the point with pixel coordinates (u, v) can be calculated.

[0041] Step 11: Based on the processed binary images at K angles described in Step 8, traverse the two-dimensional pixel coordinates of each binary image at each angle. If the current two-dimensional coordinate point (u, v) is a bubble projection point, then calculate the three-dimensional coordinates of each layer in the space to be reconstructed by back-projecting from the two-dimensional pixel coordinates (u, v) according to formula (11), that is, obtain the three-dimensional back-projection points of the current two-dimensional bubble image in the space to be reconstructed. After the traversal, take the union of all the three-dimensional back-projection points of the binary images at K angles to form a three-dimensional back-projection point set

[0042] Step 12: Based on the three-dimensional back-projection point set described in Step 11, initialize the voxel weights in the space to be reconstructed according to formula (12).

[0043]

[0044] Among them, w j is the weight factor of the jth voxel, j = 1, 2, 3,..., H, and H represents the total number of voxels in the discretized model.

[0045] Assume that the jth voxel is in the τth layer of the world coordinate system of the space to be reconstructed, then its voxel intensity x j can be initialized through formula (13).

[0046]

[0047] Among them, L τ is the number of bubble pixel points on the τ-th layer in the bubble image; is the number of back-projected points of bubbles on the τ-th layer in three-dimensional space.

[0048] Step 13: Traverse the voxels of the space to be reconstructed, judge the voxel weight of each voxel. If the weight factor w j > 0, then calculate the pixel coordinates (u, v) projected onto the image at the n-th angle from the world coordinates (X w , Y w , Z w ) of the current voxel according to formula (8) described in step 6, and accumulate the contribution values of all voxels to the point (u, v) on the projection image at the n-th angle, so as to obtain the estimated projection value t n (u, v) at the coordinate (u, v) on the projection image at the n-th angle.

[0049] Assume that the coordinate point (u, v) on the projection image is the i-th pixel on this image, then there is

[0050]

[0051] In the formula, represents the voxel weight when the projection point of the j-th voxel is the i-th pixel; i = 1, 2, 3,..., G, and G represents the total number of two-dimensional pixels.

[0052] If the weight factor w j < 0, then this point is not projected to reduce the calculation amount and improve the iteration speed.

[0053] Step 14: Calculate the mean absolute error (MAE) between the bubble images and the bubble estimated projection images at all angles using formula (15).

[0054]

[0055] Among them, I n (i) is the pixel value of the i-th pixel point of the bubble image at the n-th angle after being processed in step 8; t n (i) is the pixel value of the i-th pixel point on the bubble estimated projection image at the n-th angle obtained in step 13; i = 1, 2, 3,..., G, and G represents the total number of two-dimensional pixels.

[0056] Step 15: Judge whether to stop the iteration. If Δ e is less than error min or the number of iterations is greater than Ite m , then stop the iteration and output the result by x jThe three-dimensional voxel matrix X formed is the reconstruction result. Otherwise, go to step 16 for iterative update.

[0057] Step 16: Use the iterative formula of the MLEM method shown in formula (16) as the update criterion to perform iterative update operations.

[0058]

[0059] Where Ite is the number of iterations.

[0060] The specific operation process of the improved MLEM method based on multi-angle images is as follows:

[0061] Calculate the ratio Δp n (i) of the pixel value of the i-th pixel of the bubble image and the estimated projection image of the bubble at the n-th angle according to formula (17)

[0062] Δp n (i) = I n (i) / t n (i) Formula (17)

[0063] Then calculate the correction value V of the j-th voxel according to formula (18) j ,

[0064]

[0065] Then update the voxel intensity of the j-th voxel according to formula (19)

[0066]

[0067] After completing the update of the voxel intensity for all voxels, go to step 13 for iterative calculation. Description of the Drawings

[0068] Figure 1 : Structure diagram of multi-camera calibration

[0069] Figure 2 : Schematic diagram of bubble image acquisition

[0070] Figure 3 : Flow chart of the improved MLEM method

[0071] Figure 4 : Overall scheme diagram of bubble three-dimensional flow field reconstruction Detailed Implementation Manner

[0072] The present invention provides an improved MLEM method applied to the reconstruction of the three-dimensional flow field of bubbles. The specific implementation manner is to calibrate a multi-camera system, unify the world coordinate systems of the multi-cameras, and obtain the projection matrices from the same world coordinate system to the world coordinate systems of each camera. Among them, the calibration structure of the multi-cameras is as Figure 1 shown. Then, use multiple cameras to synchronously collect bubble images, and obtain bubble images at multiple angles at the same moment. The schematic diagram of bubble image collection is as Figure 2 shown. Finally, in view of the characteristics of backlight illumination of the three-dimensional flow field of bubbles, the method for setting the initial iteration value in the MLEM method and the projection and back-projection methods in the iterative process of the MLEM method are improved, so as to realize the reconstruction of the three-dimensional flow field of bubbles from bubble images at multiple angles according to the improved MLEM iterative method. The flow chart of the improved MLEM method is as Figure 3 shown. The overall scheme for the reconstruction of the three-dimensional flow field of bubbles is as Figure 4 shown. Its characteristics include the following steps:

[0073] The present invention designs an improved MLEM method applied to the reconstruction of the three-dimensional flow field of bubbles, and its characteristics include the following steps:

[0074] Step 1: Place K cameras around the shooting target and on the same horizontal plane, and the overall shooting angle of the cameras is less than 180 degrees.

[0075] Step 2: Place a checkerboard calibration board in the target area to be photographed, change the position of the checkerboard calibration board multiple times, and use K cameras to synchronously photograph the checkerboard calibration board to obtain multiple groups of checkerboard calibration board images at K angles at the same moment.

[0076] Step 3: Perform single calibration on each camera to obtain the internal and external parameter information of each camera.

[0077] During calibration, the computer will obtain the corner information of the checkerboard images taken by each camera, and each camera defines the coordinates of the checkerboard corner closest to the upper left corner as the origin of the world coordinates. The horizontal direction of the checkerboard is the X-axis, the vertical direction of the checkerboard is the Y-axis, and the direction perpendicular to the checkerboard plane is the Z-axis. Taking this as the world coordinate system, the external parameters will be calculated based on this world coordinate system. Since the position of the checkerboard needs to be changed during calibration and multiple images are used for calibration, the world coordinates corresponding to different images are also different, and the calculated external parameters are also different. In order to make the world coordinate systems of the K cameras have a unified reference, the external parameters of each camera are calculated based on the calibration images synchronously taken by the K cameras at a certain moment. The internal parameters of the camera represent the relationship between the camera coordinates and the pixel coordinates, and are only determined by the camera itself and can be uniquely determined.

[0078] Perform single calibration on each camera to obtain the internal parameter matrix M of the nth (1 <= n <= K) camera n , rotation matrix R n and translation vector T n . At this time, the rotation matrix and translation vector of each camera are both determined with respect to the world coordinate system of the checkerboard calibration board at the same position.

[0079] Step 4: Use camera 1 as the reference camera, and utilize the internal and external camera parameters obtained in Step 3 to calculate the position relationship parameters between other cameras and camera 1.

[0080] Perform binocular calibration on camera 1 and camera n to calculate the position relationship parameters of camera n relative to camera 1; assume that a certain point in space is q, and its world coordinate is known as q w , and the coordinates of point q in the camera coordinate systems of camera 1 and camera n are

[0081] q 1 = R 1 ·q w + T 1 Formula (1)

[0082] q n = R n ·q w + T n Formula (2)

[0083] From Formula (1) and Formula (2), the relationship between the two camera coordinates q 1 and q n is

[0084] q n = R 1n ·q 1 + T 1n Formula (3)

[0085] where

[0086]

[0087] T 1n = T n - R 1n ·T 1 Formula (5)

[0088] R 1n and T 1n are respectively the rotation matrix and translation vector between camera 1 and camera n, which are the required position relationship parameters between cameras.

[0089] Step 5: Based on the world coordinate system of Camera 1, use the position relationship matrix of Camera 1 and Camera n described in Step 4 to convert the world coordinate system of Camera n into the world coordinate system of Camera 1. Then, the external parameter matrix of Camera n is

[0090]

[0091] Step 6: The relationship between the world coordinate system (X w , Y w , Z w ) and the pixel coordinate system (u, v) is shown in Equation (7).

[0092]

[0093] Among them, Z c is the coordinate value of the current coordinate point on the Z-axis of the camera coordinate system; M is the internal parameter matrix of the camera, including f x , f y , u 0 , v 0 , s and other parameters. f x and f y represent the scale factors for the camera to transform between the pixel coordinate system and the camera coordinate system. (u 0 , v 0 ) is the corresponding coordinate of the origin of the image coordinate system in the pixel coordinate system, and s is the skew factor; N is the external parameter matrix of the camera, including the rotation matrix R and the translation vector T; MN is the camera projection matrix, and the matrix dimension is 3×4.

[0094] Substitute the internal parameters of each camera calculated in Step 3 and the external parameter matrix of each camera calculated in Step 4 into Equation (7), and the projection relationship expression between the world coordinate system and the pixel coordinate system of each camera can be obtained, as shown in Equation (8).

[0095]

[0096] In the formula, (u n , v n ) is the pixel point coordinate in the pixel coordinate system of Camera n; M n is the internal parameter matrix of Camera n.

[0097] Step 7: After completing the content described in Step 6, keep the relative positions of the cameras unchanged, place a water tank in the photographed area, introduce bubbles into the water tank, and use K cameras to synchronously photograph the bubble flow field to obtain bubble images at K angles at the same moment.

[0098] Step 8: Based on the K-angle bubble images collected in Step 7, perform denoising and binarization processing on the images to obtain the binarized bubble images I at K angles1 , I 2 ,..., I n ,..., I K .

[0099] Step 9: Set the maximum number of iterations Ite of the iterative reconstruction method m and the minimum allowable error error min .

[0100] Step 10: Represent the projection matrix MN in formula (7) described in step 6 using formula (9).

[0101]

[0102] where C ab is the matrix element in the projection matrix, a = 1, 2, 3, b = 1, 2, 3, 4.

[0103] Then the projection relationship from the world coordinate system of the camera to the pixel coordinate system can be expressed in the form shown in formula (10).

[0104]

[0105] Divide the three-dimensional space region into layers along the Z w axis of the world coordinate system. Consider the three-dimensional space region as a continuous multi-layer two-dimensional image. For the known Z w layer, its back-projection calculation formula is as shown in formula (11).

[0106]

[0107] Substitute the two-dimensional pixel coordinates (u, v) and the Z w coordinate of the current layer into the back-projection calculation formula, and the X w , Y w , Y w coordinates of the corresponding point on the Z

[0108] Step 11: Based on the processed binary images at K angles described in step 8, traverse the two-dimensional pixel coordinates of each binary image at each angle. If the current two-dimensional coordinate point (u, v) is a bubble projection point, then calculate the three-dimensional coordinates of each layer in the space to be reconstructed by back-projecting from the two-dimensional pixel coordinates (u, v) according to formula (11), that is, obtain the three-dimensional back-projection points of the current two-dimensional bubble image projected into the space to be reconstructed. After the traversal is completed, take the union of all the three-dimensional back-projection points of the binary images at K angles to form a three-dimensional back-projection point set

[0109] Step 12: Based on the three-dimensional back-projection point set described in Step 11, initialize the voxel weights of the space to be reconstructed according to Formula (12).

[0110]

[0111] where w j is the weight factor of the j-th voxel, and j = 1, 2, 3, ..., H, where H represents the total number of voxels in the discretized model.

[0112] Assume that the j-th voxel is on the τ-th layer of the world coordinate system of the space to be reconstructed. Then, its voxel intensity x can be initialized through Formula (13). j .

[0113]

[0114] where L τ is the number of bubble pixel points on the τ-th layer in the bubble image; is the number of back-projection points of bubbles on the τ-th layer in the three-dimensional space.

[0115] Step 13: Traverse the voxels of the space to be reconstructed, and judge the voxel weights of each voxel. If the weight factor w j > 0, then calculate the pixel coordinates (u, v) projected onto the image at the n-th angle from the world coordinates (X w , Y w , Z w ) of the current voxel according to Formula (8) described in Step 6, and accumulate the contribution values of all voxels to the point (u, v) on the projection image at the n-th angle, so as to obtain the estimated projection value t n (u, v) at the coordinate (u, v) on the projection image at the n-th angle.

[0116] Let the coordinate point (u, v) on the projection image be the i-th pixel of this image. Then, there is

[0117]

[0118] In the formula, represents the voxel weight when the projection point of the j-th voxel is the i-th pixel; i = 1, 2, 3, ..., G, where G represents the total number of two-dimensional pixels.

[0119] If the weight factor w j < 0, then this point is not projected to reduce the calculation amount and improve the iteration speed.

[0120] Step 14: Calculate the mean absolute error (MAE) between the bubble images and the bubble estimated projection images at all angles using Formula (15).

[0121]

[0122] Among them, I n (i) is the pixel value of the i-th pixel of the n-th angular bubble image after being processed in step 8; t n (i) is the pixel value of the i-th pixel on the estimated projection image of the bubble at the n-th angle obtained in step 13.

[0123] Step 15: Determine whether to stop the iteration. If Δ e is less than error min or the number of iterations is greater than Ite m , then stop the iteration and output the three-dimensional voxel matrix X composed of x j . X is the reconstruction result. Otherwise, go to step 16 for iterative update.

[0124] Step 16: Use the iterative formula of the MLEM method shown in formula (16) as the update criterion to perform iterative update operations.

[0125]

[0126] Among them, Ite is the number of iterations.

[0127] The specific operation process of the improved MLEM method based on multi-angle images is as follows:

[0128] Calculate the ratio Δp n (i) of the pixel values of the i-th pixel of the bubble image and the estimated projection image of the bubble at the n-th angle according to formula (17),

[0129] Δp n (i) = I n (i) / t n (i) Formula (17)

[0130] Then calculate the correction value V j of the j-th voxel by formula (18),

[0131]

[0132] Then update the voxel intensity of the j-th voxel according to formula (19),

[0133]

[0134] After completing the update of the voxel intensities for all voxels, go to step 13 for iterative calculation.

[0135] The above has schematically described the present invention and its embodiments. This description is not restrictive, and what is shown in the drawings is only one of the embodiments of the present invention. Therefore, if those of ordinary skill in the art are inspired by it and, without departing from the gist of the present invention, adopt other forms of similar components or other forms of component layout methods and creatively design technical solutions and embodiments similar to this technical solution, they shall fall within the protection scope of the present invention.

Claims

1. An improved MLEM method applied to the reconstruction of the three-dimensional bubble flow field, characterized in that, it includes the following steps: Step 1: Place K cameras around the shooting target and on the same horizontal plane, and the overall shooting angle of the cameras is less than 180 degrees; Step 2: Place a checkerboard calibration board in the target area to be photographed. Change the position of the checkerboard calibration board multiple times, and use K cameras to synchronously photograph the checkerboard calibration board to obtain multiple groups of checkerboard calibration board images at K angles at the same moment; Step 3: Perform single calibration on each camera to obtain the internal and external parameter information of each camera; During calibration, the computer will obtain the corner point information of the checkerboard image captured by each camera. Each camera defines the coordinates of the checkerboard corner point closest to the upper left corner as the origin of the world coordinates. The horizontal direction of the checkerboard is the X-axis, the vertical direction of the checkerboard is the Y-axis, and the direction perpendicular to the checkerboard plane is the Z-axis. Taking this as the world coordinate system, the external parameters will be calculated based on this world coordinate system; Since the position of the checkerboard needs to be changed during the calibration process and multiple images are used for calibration, the world coordinates corresponding to different images are also different, and the calculated external parameters are also different; In order to make the world coordinate systems of K cameras have a unified reference, the external parameters of each camera are calculated based on the calibration images synchronously captured by K cameras at a certain moment; The internal parameters of the camera represent the relationship between the camera coordinates and the pixel coordinates and are only determined by the camera itself and can be uniquely determined; Perform single calibration on each camera to obtain the internal parameter matrix M of the nth camera n , rotation matrix R n and translation vector T n . At this time, the rotation matrix and translation vector of each camera are both determined with respect to the world coordinate system of the checkerboard calibration board at the same position; Step 4: Use camera 1 as the reference camera, and use the internal and external parameters of the camera obtained in Step 3 to calculate the position relationship parameters between other cameras and camera 1; Perform binocular calibration on Camera 1 and Camera n, and calculate the position relationship parameters of Camera n relative to Camera 1; assume that a certain point in space is q, and its world coordinates are known as q w , and the coordinates of point q in the camera coordinate systems of Camera 1 and Camera n are q 1 = R 1 ·q w + T 1 Equation (1) q n = R n ·q w + T n Equation (2) From formulas (1) and (2), the relationship between the two camera coordinates q 1 and q n is q n = R 1n ·q 1 + T 1n Formula (3) Among them, T 1n = T n - R 1n · T 1 Equation (5) R 1n and T 1n are respectively the rotation matrix and the translation vector between Camera 1 and Camera n, which are the parameters of the inter-camera position relationship to be obtained; Step 5: Based on the world coordinate system of camera 1, use the position relationship matrix between camera 1 and camera n described in Step 4 to convert the world coordinate system of camera n to the world coordinate system of camera 1, then the external parameter matrix of camera n is Step 6: The relationship between the world coordinate system (X w , Y w , Z w ) and the pixel coordinate system (u, v) is shown in Equation (7); Among them, Z c is the coordinate value of the current coordinate point on the Z-axis of the camera coordinate system; M is the internal parameter matrix of the camera, including f x , f y , u 0 , v 0 , s, f x and f y represent the scale factors for the conversion between the pixel coordinate system and the camera coordinate system of the camera. (u 0 , v 0 ) is the corresponding coordinate of the origin of the image coordinate system in the pixel coordinate system, and s is the skew factor; N is the external parameter matrix of the camera, including the rotation matrix R and the translation vector T; MN is the camera projection matrix, and the matrix dimension is 3×4; Substitute the internal parameters of each camera calculated in Step 3 and the external parameter matrices of each camera calculated in Step 4 into formula (7), and the projection relationship expression between the world coordinate system and the pixel coordinate systems of each camera can be obtained, as shown in formula (8); where (u n , v n ) are the pixel coordinates in the pixel coordinate system of camera n; M n is the internal parameter matrix of camera n; Step 7: After completing Step 6, keep the relative positions of the cameras unchanged, place a water tank in the photographed area, introduce bubbles into the water tank, and use K cameras to synchronously photograph the bubble flow field to obtain bubble images at K angles at the same moment; Step 8: Based on the K angular bubble images collected in Step 7, perform denoising and binarization processing on the images to obtain K binarized bubble images I 1 , I 2 , …, I n , …, I K ; Step 9: Set the maximum number of iterations Ite of the iterative reconstruction method m and the minimum allowable error error min ; Step 10: Represent the projection matrix MN in formula (7) described in Step 6 with formula (9); Among them, C ab is the matrix element in the projection matrix, where a = 1, 2, 3, b = 1, 2, 3, 4; Then the projection relationship from the world coordinate system of the camera to the pixel coordinate system can be expressed in the form shown in formula (10); The three-dimensional space region is stratified along the Z-axis of the world coordinate system, and the three-dimensional space region is regarded as a continuous multi-layer two-dimensional image. For the known Z layer, its back-projection calculation formula is shown in formula (11); w axis, regarding the three-dimensional space region as a continuous multi-layer two-dimensional image, for the known Z w layer, its back-projection calculation formula is as shown in formula (11); Substitute the two-dimensional pixel coordinates (u, v) and the Z coordinate of the current layer w into the back-projection calculation formula, and the X w , Y w coordinates of the corresponding point on the Z w layer of the point with pixel coordinates (u, v) can be calculated; Step 11: Based on the binarized images at K angles processed in Step 8, traverse the two-dimensional pixel coordinates of each binarized image at an angle. If the current two-dimensional coordinate point (u, v) is a bubble projection point, calculate the three-dimensional coordinates of each layer in the space to be reconstructed back-projected from the two-dimensional pixel coordinates (u, v) according to formula (11), that is, obtain the three-dimensional back-projection points of the current two-dimensional bubble image back-projected into the space to be reconstructed. After the traversal, take the union of all the three-dimensional back-projection points of the binarized images at K angles to form a three-dimensional back-projection point set Step 12: Based on the three-dimensional back-projection point set described in Step 11, initialize the voxel weights of the space to be reconstructed according to formula (12); where, w j is the weight factor of the j-th voxel, j = 1, 2, 3, ..., H, and H represents the total number of voxels in the discretized model; Assume that the j-th voxel is in the τ-th layer of the world coordinate system of the space to be reconstructed, then its voxel intensity x can be initialized by formula (13). j ; where L τ is the number of bubble pixel points on the τ-th layer in the bubble image; is the number of back-projected points of the bubbles on the τ-th layer in the three-dimensional space; Step 13: Traverse the voxels of the space to be reconstructed, and judge the voxel weights of each voxel. If the weight factor w j > 0, calculate the pixel coordinates (u, v) projected onto the image of the nth angle from the world coordinates (X w , Y w , Z w ) of the current voxel according to formula (8) described in step 6, and accumulate the contribution values of all voxels to the point (u, v) on the projection image of the nth angle, so as to obtain the estimated projection value t n (u, v) at the coordinate (u, v) on the projection image of the nth angle; Let the coordinate point (u, v) on the projection image be the i-th pixel on this image, then there is In the formula, represents the voxel weight when the projection point of the j-th voxel is the i-th pixel; i = 1, 2, 3, ..., G, where G represents the total number of two-dimensional pixels; If the weight factor w j < 0, then this point is not projected to reduce the computational amount and improve the iteration speed; Step 14: Calculate the mean absolute error (MAE) of the bubble images and the bubble estimated projection images at all angles using formula (15); Among them, I n (i) is the pixel value of the i-th pixel point of the n-th angular bubble image after the processing in step 8; t n (i) is the pixel value of the i-th pixel point on the estimated projection image of the bubble at the n-th angle obtained in step 13; i = 1, 2, 3,..., G, where G represents the total number of two-dimensional pixels; Step 15: Determine whether to stop iteration; if Δ e is less than error min or the number of iterations is greater than Ite m , then stop iteration and output the three-dimensional voxel matrix X composed of x j , and X is the reconstruction result; otherwise, go to Step 16 for iterative update; Step 16: Use the MLEM method iteration formula shown in formula (16) as the update criterion to perform iterative update operations; where Ite is the number of iterations; The specific operation process of the improved MLEM method based on multi-angle images is as follows: Calculate the ratio Δp of the pixel value of the i-th pixel of the bubble image and the estimated projection image of the bubble at the n-th angle according to formula (17). n (i). Δp n (i) = I n (i) / t n (i) Equation (17) Then, the correction value V of the j-th voxel is calculated by formula (18). j , Then, update the voxel intensity of the j-th voxel according to formula (19), After completing the update of the voxel intensity for all voxels, go to step 13 for iterative calculation.

Images