Convolution reconstruction method, system, equipment and medium for tomographic particle image velocimetry

Abstract

The present invention relates to a convolution reconstruction method, system, device, and medium for tomographic particle image velocimetry. The convolution reconstruction method for tomographic particle image velocimetry comprises: S101 obtaining a particle image required for reconstruction and setting reconstruction parameters; S102 initializing the discretized spatial voxels; S103 re-dividing the original pixel grid of the particle image into levels to obtain a high-resolution pixel grid; S104 obtaining the grayscale of the high-resolution particle image; S105 calculating an approximate weight matrix of the spatial voxels based on the division levels; S106 projecting the sphere centers of all non-zero voxels onto a high-resolution pixel plane according to a mapping function, and assigning a value to the nearest pixel center to obtain a particle projection matrix; S107 using the approximate weight matrix as a convolution kernel and performing convolution with the particle projection matrix to generate a projection image; S108 calculating the ratio of each pixel corresponding to the projection image and the high-resolution image; and S109 updating the grayscale of the spatial voxels. The present invention improves the efficiency and quality of particle field reconstruction for tomographic particle image velocimetry.

Description

Technical Field

[0001] The present invention relates to the field of experimental fluid mechanics technology, and in particular to a convolution reconstruction method, system, equipment and medium for tomographic particle image velocimetry. Background Art

[0002] Tomographic Particle Image Velocimetry (Tomo-PIV) is a key method for quantitatively measuring three-dimensional flow fields in experimental fluid dynamics. This technique records particle motion in the flow field from multiple perspectives, uses advanced reconstruction algorithms to construct the three-dimensional particle distribution of the flow field, and finally calculates the three-dimensional velocity field distribution using a three-dimensional cross-correlation algorithm. The advent of Tomo-PIV not only marks a significant advancement in flow measurement technology from two-dimensional to three-dimensional, but also provides a more powerful tool for studying complex flow phenomena such as turbulence and eddies.

[0003] The core step of Tomo-PIV technology is spatial grayscale reconstruction (also known as particle field reconstruction). Spatial grayscale reconstruction involves finding the spatial grayscale distribution E given a known particle image P and a weight function W. This is the inverse process of projection imaging, where WE = P. Because the number of discrete voxels far exceeds the known number of pixels, this inverse problem has an indeterminate solution. Currently, the most common reconstruction algorithm is the multiplicative algebraic reconstruction technique (MART), whose iterative formula is:

[0004]

[0005] Where (x i ,y i ) is the coordinate of pixel i, P(x i ,y i ) is the pixel grayscale of pixel i, (X j , Y j , Z j ) is the coordinate of voxel j, B(X j , Y j , Z j ) is the voxel grayscale of voxel j, w i，j is the weight coefficient, representing the contribution of voxel j to the grayscale of pixel i. This contribution is related to the projection distance; smaller distances result in greater contributions. The weight coefficient ranges from 0 to 1. k represents the number of iterations, and μ is the relaxation coefficient.

[0006] While the MART method can effectively reconstruct particle fields in three dimensions, its greatest limitation lies in its low computational efficiency. The main reasons for this are: first, the MART method calculates the weight coefficient for each voxel individually, and the computational effort is proportional to the cube of the reconstructed volume size; second, the MART method requires a large amount of memory and places high demands on the computing equipment. These issues have severely hampered the practical application of Tomo-PIV.

[0007] In response to the limitations of the MART method, scholars have conducted continuous research and proposed a variety of improved fast reconstruction algorithms. Among them, multiplicative line of sight-simultaneous MART (MLOS-SMART) can significantly reduce memory usage while improving reconstruction speed, making it the most widely used fast reconstruction algorithm. However, this algorithm still does not achieve the fastest reconstruction speed and still takes a long time to process large-scale data sets. In most cases, it can only quantitatively calculate the instantaneous three-dimensional flow field and cannot calculate the time-resolved flow field in large quantities. How to improve the efficiency of Tomo-PIV's particle field reconstruction is an urgent problem that needs to be solved. Summary of the Invention

[0008] The present invention provides a convolution reconstruction method, system, device and medium for tomographic particle image velocimetry to solve the problem of low efficiency of particle field reconstruction in Tomo-PIV.

[0009] To achieve the above objectives, in a first aspect, the present invention provides a convolution reconstruction method for tomographic particle image velocimetry, comprising:

[0010] S101 obtains particle images required for reconstruction and sets reconstruction parameters, wherein the particle images are multiple particle images obtained by taking multiple cameras set at different angles;

[0011] S102 initializes the discretized spatial voxels;

[0012] S103: re-dividing the original pixel grid of the particle image based on the division level to obtain a high-resolution pixel grid, wherein the original pixel size in the high-resolution pixel grid is converted into an interpolated pixel size;

[0013] S104 interpolates the particle image on the original pixel plane to the high-resolution pixel grid to obtain a high-resolution particle image grayscale image.

[0014] S105 Calculate the approximate weight matrix of the spatial voxels according to the division level

[0015] S106 projects the sphere centers of all non-zero voxels onto the high-resolution pixel plane according to the mapping function, finds the pixel centers closest to them, assigns the initial grayscale values ​​of the non-zero voxels to the corresponding pixels, and obtains a particle projection matrix;

[0016] S107 uses the approximate weight matrix as a convolution kernel and performs convolution with the particle projection matrix to generate a projection image

[0017] S108 Exploitation Calculate the projection image With the high-resolution image The ratio of each pixel

[0018] S109 Utilization Formula Update the grayscale of the spatial voxel, and the grayscale of the spatial voxel is expressed by E(X j , Y j , Z j ) means that E(X j , Y j , Z j ) is the coordinate of the voxel in the j-space volume.

[0019] Preferably, the multiple cameras are four cameras distributed orthogonally in an X shape, the angle between two adjacent cameras is 60 degrees, and the reconstruction parameters include at least a mapping function coefficient, an iteration number, and a threshold coefficient;

[0020] The method further includes: S110 repeatedly executing all steps S106 to S109 until a preset number of iterations of the particle image is reached.

[0021] Preferably, the division level is 4.

[0022] Preferably, the approximate weight matrix of the voxels is calculated according to the division level Specifically, assuming that each spatial voxel is projected to the pixel center, an approximate weight matrix is ​​calculated, specifically:

[0023] The voxel simplified into a sphere is projected onto the high-resolution pixel plane, and the circle is called the voxel projection circle;

[0024] Calculate the area A of the high-resolution pixels overlapping the projected circle mn , m, z represent the indices of the overlapping pixels. If the pixels are completely within the projection circle, the overlapping area is d pixel 2 / R 2 If there is partial overlap, it can be approximated by the Monte Carlo method, where d pixel is the original pixel size, R represents the division level, and the optimal value of R is 4;

[0025] Calculate the approximate weight matrix The calculation formula is A mn / A projection , where A projectoin The area of ​​the voxel projection circle is represented by the formula Calculated;

[0026] According to the mapping function, the sphere centers of all non-zero voxels are projected onto the high-resolution pixel plane, the pixel centers closest to them are found, and the initial grayscale values ​​of these non-zero voxels are assigned to the corresponding pixels to obtain the particle projection matrix. If there are multiple voxels corresponding to one pixel, the grayscale values ​​of these voxels are accumulated and assigned to the pixel.

[0027] Preferably, the formula Update the grayscale of the spatial voxel. The steps are:

[0028] Residual Find the logarithm and get

[0029] Using the approximate weight matrix right Perform convolution operation to obtain

[0030] Exponentiate the result of the convolution operation to obtain Where r represents the pixel that overlaps with the projected circle of voxel j;

[0031] Will and Multiply and update the grayscale of the spatial voxel.

[0032] Preferably, the mapping function is:

[0033] F(x)=a0+a1X+a2Y+a3X 2 +a4XY+a5Y 2 +a6X 3 +a7X2Y+a8XY 2 +a9Y 3

[0034]

[0035] To achieve the above objectives, the second aspect is for reconstructing a three-dimensional flow field, including:

[0036] A reconstruction setting module is used to obtain particle images required for reconstruction and set reconstruction parameters, wherein the particle images are multiple particle images obtained by taking multiple cameras set at different angles;

[0037] Discrete initialization module, used to initialize the discretized spatial voxels;

[0038] a high-resolution pixel grid acquisition module, configured to re-divide the original pixel grid of the particle image based on the division levels to obtain a high-resolution pixel grid, wherein the original pixel size in the high-resolution pixel grid is converted into an interpolated pixel size;

[0039] The interpolation module is used to interpolate the particle image on the original pixel plane to the high-resolution pixel grid to obtain the high-resolution particle image grayscale

[0040] An approximate weight matrix calculation module is used to calculate the approximate weight matrix of the spatial voxels according to the division level

[0041] The particle projection matrix acquisition module is used to project the sphere centers of all non-zero voxels onto the high-resolution pixel plane according to the mapping function, find the pixel centers closest to them, assign the initial grayscale values ​​of the non-zero voxels to the corresponding pixels, and obtain the particle projection matrix.

[0042] A projection image generation module is used to use the approximate weight matrix As the convolution kernel and particle projection matrix Perform convolution to generate the projected image

[0043] Pixel weight calculation module, used to use Calculate the projection image With the high-resolution particle image grayscale The ratio of each pixel

[0044] Update the spatial voxel grayscale module to use the formula Update the grayscale of the spatial voxel. The grayscale of the voxel of the spatial particle field is updated with E(X j , Y j , Z j ) means, (X j , Y j , Z j ) is the coordinate of voxel j.

[0045] Preferably, the multiple cameras are four cameras distributed orthogonally in an X shape, and the angle between two adjacent cameras is 60 degrees. The reconstruction parameters include at least a mapping function coefficient, an iteration number, and a threshold coefficient; and further include: an iteration module for repeatedly executing the particle projection matrix acquisition module, the projection image generation module, the pixel proportion calculation module, and the update space voxel grayscale module until the preset iteration number of the particle image is reached.

[0046] Preferably, the approximate weight matrix calculation module is specifically used to:

[0047] Assuming that each spatial voxel is projected to the pixel center, the approximate weight matrix is ​​calculated once:

[0048] The voxel simplified into a sphere is projected onto the high-resolution pixel plane, and the circle is called the voxel projection circle;

[0049] Calculate the area A of the high-resolution pixels overlapping the projected circle mn , m, n represent the indices of the overlapping pixels. If the pixels are completely within the projection circle, the overlapping area is d pixel 2 / R 2 ; If there is partial overlap, it can be approximated by the Monte Carlo method;

[0050] Calculate the approximate weight matrix The calculation formula is A mn / A projection , where A projection The area of ​​the voxel projection circle is represented by the formula Calculated;

[0051] According to the mapping function, the sphere centers of all non-zero voxels are projected onto the high-resolution pixel plane, the pixel centers closest to them are found, and the initial grayscale values ​​of these non-zero voxels are assigned to the corresponding pixels to obtain the particle projection matrix. If there are multiple spatial voxels corresponding to one pixel, the grayscale values ​​of these spatial voxels are accumulated and assigned to the pixel.

[0052] Preferably, the updating spatial voxel grayscale module is specifically used to:

[0053] Residual Find the logarithm and get

[0054] Using the approximate weight matrix right Perform convolution operation to obtain

[0055] Exponentiate the result of the convolution operation to obtain Where r represents the pixel that overlaps with the projected circle of voxel j;

[0056] Will and Multiply and update the grayscale of the spatial voxel.

[0057] To achieve the above-mentioned purpose, in a third aspect, the present invention also relates to an electronic device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that when the processor executes the computer program, the convolution reconstruction method for tomographic particle image velocimetry mentioned above is implemented.

[0058] To achieve the above objectives, in a fourth aspect, the present invention further relates to a computer-readable storage medium, in which instructions are stored, and when the instructions are run, the above-mentioned convolution reconstruction method for tomographic particle image velocimetry is executed.

[0059] The present invention relates to a convolution reconstruction method, system, device and medium for tomographic particle image velocimetry, which has the following beneficial effects compared with the prior art:

[0060] Based on the MLOS-SMART method, a convolutional multiplication algebraic reconstruction algorithm (Convolutional MART, or Conv-MART) was proposed to improve the efficiency of particle field reconstruction using Tomo-PIV. This algorithm re-divides the pixel grid to obtain high-resolution pixels and replaces the weight matrices with approximate weight matrices to reduce computation time. Furthermore, the convolution operation is introduced based on the SMART method to further accelerate the computation process, resulting in the Conv-MART algorithm. Compared with the MLOS-SMART method, this algorithm has the following advantages: 1) It reduces the computation time of the weight matrix and reduces the number of computational steps, significantly improving reconstruction efficiency; 2) Conv-MART achieves higher reconstruction quality than MLOS-SMART when the particle image contains noise or the particle concentration is high. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 This is a flow chart of a convolution reconstruction method for tomographic particle image velocimetry in Example 1 of the present invention;

[0062] Figure 2 Schematic diagram of pixel grid division for a convolution reconstruction method for tomographic particle image velocimetry in Example 1 of the present invention;

[0063] Figure 3 This is a flow chart of step 105 of a convolution reconstruction method for tomographic particle image velocimetry in Example 1 of the present invention;

[0064] Figure 4This is a flow chart of step 109 of a convolution reconstruction method for tomographic particle image velocimetry in Example 1 of the present invention;

[0065] Figure 5 This is a reconstruction time comparison diagram of Example 1 of a convolution reconstruction method for tomographic particle image velocimetry in Example 1 of the present invention;

[0066] Figure 6 This is a velocity field diagram of Example 1 of a convolution reconstruction method for tomographic particle image velocimetry in Example 1 of the present invention;

[0067] Figure 7 This is a reconstruction time comparison diagram of Example 2 of a convolution reconstruction method for tomographic particle image velocimetry in Example 1 of the present invention;

[0068] Figure 8 This is a velocity field diagram of Example 2 of a convolution reconstruction method for tomographic particle image velocimetry in Example 1 of the present invention;

[0069] Figure 9 Schematic diagram of the structure of a convolution reconstruction system for tomographic particle image velocimetry in the second embodiment of the present invention;

[0070] Figure 10 It is a structural diagram of an electronic device in embodiment 3 of the present invention. DETAILED DESCRIPTION

[0071] The present invention will be further described in detail below with reference to the accompanying drawings and examples. It will be understood that the specific embodiments described herein are intended only to illustrate the present invention and are not intended to limit the present invention. It should also be noted that, for ease of description, the accompanying drawings only illustrate portions relevant to the present invention, not all structures.

[0072] Example 1

[0073] A convolution reconstruction method for tomographic particle image velocimetry, see Figures 1-8 , for the reconstruction of three-dimensional flow field, the present invention proposes a convolution multiplication algebraic reconstruction method (Conv-MART) to improve the efficiency of particle field reconstruction of Tomo-PIV. In the present invention, the pixel grayscale of the particle image is expressed as P(x i ,y i ) means, (x i ,y i ) is the coordinate of the pixel. The voxel grayscale of the spatial particle field is expressed by E(X j , Y j , Z j ) means, (X j , Y j , Z j ) is the coordinate of voxel j. The weight coefficient between pixels and voxels is w i，j. Voxels are simplified to spheres, and pixels remain squares. The present invention mainly redivides the particle image into a high-resolution pixel grid, and calculates an approximate weight matrix to replace different weight matrices. On this basis, convolution operations are used for iterative calculations to obtain a reconstructed particle field. In addition, the present invention can easily utilize an image processing unit (GPU) to further accelerate the reconstruction.

[0074] The process includes the following steps: S101 to S110.

[0075] S101: Obtain particle images required for reconstruction and set reconstruction parameters, wherein the particle images are multiple particle images obtained by taking multiple cameras set at different angles.

[0076] In this embodiment, the multiple cameras are four cameras arranged in an X-shaped orthogonal pattern, with the angle between two adjacent cameras being 60 degrees. The reconstruction parameters include at least the mapping function coefficients, the number of iterations, and the threshold coefficient. If a GPU is required, the particle images and mapping function coefficients in this step must be initialized to GPU data.

[0077] S102: Initialize the discretized spatial voxels.

[0078] The initialization method includes but is not limited to uniform initialization, MLOS initialization or other initialization methods. If the MLOS initialization method is used, the spatial voxels need to be assigned values ​​according to the product of the pixel grayscales of different cameras.

[0079] S103: Re-dividing the original pixel grid of the particle image based on the division level to obtain a high-resolution pixel grid, wherein the original pixel size in the high-resolution pixel grid becomes an interpolated pixel size.

[0080] The pixel size is determined by the original d pixel (usually 1 pixel) becomes the divided d pixel / R, R represents the resize level. Figure 2 The pixel grid is divided into two parts. The pixel size after division is called interpolation pixel size. The present invention has been tested and found that the value of R is 4, which is the most suitable. At this time, the interpolation pixel size is d pixel / 4.

[0081] S104 interpolates the particle image on the original pixel plane to the high-resolution pixel grid to obtain a high-resolution particle image grayscale image.

[0082] Among them, interpolation methods include but are not limited to bilinear interpolation, polynomial interpolation, adjacent point interpolation, etc.

[0083] S105: Calculate the approximate weight matrix of the spatial voxels according to the division level

[0084] In this step, it is assumed that each voxel in space is projected to the pixel center. No projection process is required. The approximate weight matrix calculated for different voxels is the same, so it only needs to be calculated once. The higher the resolution of the generated image, the more this assumption is true. Figure 3 This is the flow chart of this step, and the specific process includes S1051-S1053:

[0085] Assuming that each spatial voxel is projected to the pixel center, the approximate weight matrix is ​​calculated once, specifically:

[0086] S1051: Projecting the voxel simplified into a sphere onto a high-resolution pixel plane. This circle is called a voxel projection circle.

[0087] S1052: Calculate the area A of the high-resolution pixel overlapping with the projected circle mn , m, n represent the indices of the overlapping pixels. If the pixels are completely within the projection circle, the overlapping area is d pixel 2 / R 2 If there is partial overlap, it can be approximated by the Monte Carlo method, where d pixel is the original pixel size, R represents the division level, and the optimal value of R is 4.

[0088] S1053: Calculate approximate weight matrix The calculation formula is A mn / A projection , where A projection The area of ​​the voxel projection circle is represented by the formula Calculated.

[0089] matrix For square matrices, element values ​​and matrix dimensions are related to the partitioning levels R and d voxel with d pixel It is related to the ratio M.

[0090] S106: Project the sphere centers of all non-zero voxels onto the high-resolution pixel plane according to the mapping function, find the pixel centers closest to them, assign the initial grayscale values ​​of these non-zero voxels to the corresponding pixels, and obtain the particle projection matrix If there are multiple voxels corresponding to one pixel, the grayscale values ​​of these voxels are accumulated and assigned to the pixel.

[0091] In this embodiment, the mapping function used is formula (1, 2), and a pinhole imaging model or other relationship formulas may also be used.

[0092] F(x)=a0+a1X+a2Y+a3X 2 +a4XY+a5Y 2 +a6X 3 +a7X 2 Y+a8XY 2 +a9Y 3 (1)

[0093] F(y)=b0+b1X+b2Y+b3X 2 +b4XY+b5Y 2 +b6X 3 +b7X 2 Y+b8XY 2 +b9Y 3 (2).

[0094] S107: Using the approximate weight matrix as a convolution kernel and the particle projection matrix to perform convolution to generate a projection image

[0095] Specifically: Use the convolution formula to generate the projection image The calculation formula for this step is as follows:

[0096]

[0097] Among them, "*" indicates that the matrix Perform convolution operation, and the convolution kernel is the weight matrix obtained in step 105 Generate projection image using convolution formula The calculation formula for this step is as follows:

[0098]

[0099] Among them, "*" indicates that the matrix Perform convolution operation, and the convolution kernel is the weight matrix obtained by S105

[0100] S108: Utilization Calculate the projection image With the high-resolution image The ratio of each pixel

[0101] S109: Update the grayscale of the spatial voxel using formula (3), and the grayscale of the spatial voxel is expressed by E(X j , Y j , Z j ) means, (X j , Y j , Z j ) is the coordinate of the voxel in the j-space volume.

[0102]

[0103] Formula (3) is equivalent to the SMART algorithm, but has higher computational efficiency.

[0104] In this embodiment, S109 specifically includes S1091-S1094:

[0105] S1091: Residual Find the logarithm and get

[0106] S1092: Using the approximate weight matrix right Perform convolution operation to obtain

[0107] S1093: Exponentiate the result obtained by the convolution operation to obtain where r represents the pixels that overlap with the projected circle of voxel j.

[0108] It can be mathematically proven Equivalent to the SMART method And it is more computationally efficient. Where r represents the pixels that overlap with the projected circle of voxel j.

[0109] S1094: and Multiply and update the grayscale of the spatial voxel.

[0110] S110: Repeat all the above steps S106-S109 until the preset number of iterations of the particle image is reached.

[0111] The above reconstruction process does not involve complex logical judgments, only simple algebraic operations, and can be fully implemented on the GPU, thereby further improving the reconstruction efficiency. Through the above steps, the reconstructed particle field is obtained in a very short time while ensuring the reconstruction quality.

[0112] In order to better illustrate the invention, two examples are given below:

[0113] Example 1:

[0114] Example 1 tests the performance of the present invention by reconstructing the particle field generated by the simulation. In this example, we simulated a vortex ring with a voxel number of 700×700×140voxels and a resolution of 0.05mm / voxel. 24,500 tracer particles were randomly distributed, and the grayscale of each spatial particle obeyed the Gaussian distribution. Then we projected the particles in the vortex ring onto four cameras to obtain four particle images. The four cameras are distributed in an X shape, and the angle between two adjacent cameras is 60°. The resolution of the particle image is 900×900 pixels, and the pixel size is 0.05mm. We used the Conv-MART method and the MLOS-SMART method to reconstruct the particle field at two consecutive moments in the vortex ring, and tested the performance of the present invention by comparing the two methods. The specific steps are as follows:

[0115] Step 1: Randomly generate a particle field with a concentration of 0.05ppp (particles per pixel) within a 700×700×140 voxel area. This particle field is the particle field at time A. Then, calculate the particle displacement after a short period of time based on the velocity of the vortex ring to obtain the particle field at time B.

[0116] Step 2: Input the particle field at time A and project it onto four orthogonally distributed camera planes, with the angle between the sight lines of adjacent cameras being 60 degrees.

[0117] Step 3: Use the MLOS method to initialize the grayscale of the particle field and find all non-zero voxels.

[0118] Step 4: Set the division level R to 4, the interpolation pixel size to 0.05 / 4mm, and re-divide the pixel plane to obtain a high-resolution pixel plane.

[0119] Step 5: Use bilinear interpolation to interpolate the particle image on the high-resolution pixel plane to obtain

[0120] Step 6: Project the voxels simplified into a sphere onto the high-resolution pixels. Use the Monte Carlo method to calculate the overlapping area A between the high-resolution pixels and the projected circle. mn , and then calculate the ratio of the overlapping area to the projected circle area to obtain the approximate weight matrix

[0121] Step 7: Use mapping functions (1) and (2) to calculate the particle projection matrix Reusing the convolution formula Generate projection image

[0122] Step 8: Use the formula To calculate the projected image With high-resolution images The ratio of each pixel between .

[0123] Step 9: Use formula (3) to reconstruct the particle field and set the number of iterations to 5.

[0124] Step 10: Input the particle field at time B, keep the other reconstruction parameters unchanged, repeat steps 2 to 9, and reconstruct the particle field at time B.

[0125] Step 11: Use the GPU version of the algorithm and repeat steps 2 to 10 to reconstruct the particle field. This step is used to test the acceleration performance of the Conv-MART method combined with GPU. The final calculation time is compared with the results of Figure 5 As shown in Figure 2, the average time for reconstructing a particle field by the MLOS-SMART method, the Conv-MART method, and the GPU-Conv-MART method is shown respectively.

[0126] Step 12: The particle field reconstructed by the Conv-MART method is calculated using the multi-layer iterative volume deformation method to obtain the velocity field. The result is as follows: Figure 6 shown.

[0127] The only difference between the Conv-MART method and the GPU-Conv-MART method is the reconstruction time. The reconstructed particle field and the calculated velocity field are the same. Figure 6 Only the velocity field obtained by the Conv-MART method is shown. The arrows in the figure represent the velocity vectors obtained by the solution, and the purple isosurfaces are the Q-criterion vorticity isosurfaces.

[0128] Example 2:

[0129] The second example involves a left ventricular Tomo-PIV experiment, primarily used to measure blood flow within the left ventricle after the native mitral valve has been replaced with a mechanical valve. In this experiment, four high-speed cameras equipped with 100mm / f2.8 lenses were arranged in a line, sampling at 200Hz, with a 36° angle between adjacent cameras. The cameras had a resolution of 2016×2016, with a pixel size of 0.05mm. The tracer particles used in the experiment were 50μm-diameter rhodamine B fluorescent particles with a particle concentration of approximately 0.045ppb. Furthermore, to simulate left ventricular blood flow, the solution used was a water-glycerol mixture with a volume ratio of 62% water to 38% glycerol. This solution had a refractive index of 1.395 and a viscosity of 3.86mPa·s at 25.6°C. The reconstructed measurement volume is 65 × 70 × 57 mm³, containing 1300 × 1400 × 1100 voxels. In this example, we use these four high-speed cameras to capture continuous particle images in a time series. The particle images at the first two moments (defined as A and B) are reconstructed using the MLOS-SMART and Conv-MART methods.

[0130] Step 1: Input the particle image at time A, initialize it using the MLOS method, and find all non-zero voxels.

[0131] Step 2: Use the MLOS method to initialize the grayscale of the particle field and find all non-zero voxels.

[0132] Step 3: Set the division level R to 5, the interpolation pixel size to 0.01 mm, and re-divide the pixel plane to obtain a high-resolution pixel plane.

[0133] Step 4: Use bilinear interpolation to interpolate the particle image on the high-resolution pixel plane to obtain

[0134] Step 5: Project the voxels simplified into a sphere onto the high-resolution pixels, and overlap the center of the projected circle with the center of the nearest pixel. Use the Monte Carlo method to calculate the overlapping area A between the pixel and the projected circle. mn , and then calculate the ratio of the overlapping area to the projected circle area to obtain the approximate weight matrix

[0135] Step 6: Use mapping functions (1) and (2) to calculate the particle projection matrix Reusing the convolution formula Generate projection image

[0136] Step 7: Use the formula To calculate the projected image With high-resolution images The ratio of each pixel between .

[0137] Step 8: Use formula (3) to reconstruct the particle field and set the number of iterations to 5.

[0138] Step 9: Input the particle image at time B, keep the other reconstruction parameters unchanged, repeat steps 1 to 8, and reconstruct the particle field at time B.

[0139] Step 10: Using the GPU version of the algorithm, repeat steps 1 to 9 to reconstruct the particle field at time AB. This step is used to test the acceleration performance of the Conv-MART method combined with GPU. The final calculation time is compared with the results of Figure 7 As shown in Figure 2, the average time for reconstructing a particle field by the MLOS-SMART method, the Conv-MART method, and the GPU-Conv-MART method is shown respectively.

[0140] Step 11: The particle field reconstructed by the Conv-MART method is calculated using the multi-layer iterative volume deformation method to obtain the velocity field. The result is as follows: Figure 8 shown.

[0141] The only difference between the Conv-MART method and the GPU-Conv-MART method is the reconstruction time. The reconstructed particle field and the calculated velocity field are the same. Figure 8 Only the velocity field obtained by the Conv-MART method is shown. The arrows in the figure represent the velocity vectors obtained by the solution, and the green isosurfaces are the Q-criterion vorticity isosurfaces.

[0142] S10: Example 2

[0143] A convolution reconstruction system for tomographic particle image velocimetry, used for Figures 1-8, for the reconstruction of three-dimensional flow field, the present invention proposes a convolution multiplication algebraic reconstruction method (Conv-MART) to improve the particle field reconstruction efficiency of Tomo-PIV. In the present embodiment, voxels are simplified to spheres and pixels are still squares. This system mainly obtains a high-resolution pixel grid by re-dividing the particle image, and calculates an approximate weight matrix to replace different weight matrices. On this basis, convolution operation is used to perform iterative calculation to obtain the reconstructed particle field. In addition, this system can easily utilize an image processing unit (GPU) to further accelerate the reconstruction. Among them, this system is implemented using an electronic device with a central processing unit, such as a PC, server, intelligent terminal, GPU cluster, etc. This system includes a reconstruction setting module 71, a discrete initialization module 72, a high-resolution pixel grid acquisition module 73, an interpolation module 74, an approximate weight matrix calculation module 75, a particle projection matrix acquisition module 76, a projection image generation module 77, a pixel proportion calculation module 78, an update space voxel grayscale module 79 and an iteration module 80.

[0144] The reconstruction setting module 71 is used to obtain particle images required for reconstruction and set reconstruction parameters, wherein the particle images are multiple particle images obtained by taking pictures with multiple cameras set at different angles.

[0145] The multiple cameras are four cameras distributed orthogonally in an X shape, and the angle between two adjacent cameras is 60 degrees. The reconstruction parameters at least include a mapping function coefficient, an iteration number, and a threshold coefficient.

[0146] The discrete initialization module 72 is used to initialize the discretized spatial voxels.

[0147] The high-resolution pixel grid acquisition module 73 is used to re-divide the original pixel grid of the particle image based on the division level to obtain a high-resolution pixel grid, in which the original pixel size in the high-resolution pixel grid becomes an interpolated pixel size.

[0148] The interpolation module 74 is used to interpolate the particle image on the original pixel plane to the high-resolution pixel grid to obtain the high-resolution particle image grayscale.

[0149] The approximate weight matrix calculation module 75 is used to calculate the approximate weight matrix of the spatial voxels according to the division level.

[0150] The particle projection matrix acquisition module 76 is used to project the sphere centers of all non-zero voxels onto the high-resolution pixel plane according to the mapping function, find the pixel centers closest to them, assign the initial grayscale values ​​of the non-zero voxels to the corresponding pixels, and obtain the particle projection matrix

[0151] Projection image generation module 77, used to use the approximate weight matrix As the convolution kernel and particle projection matrix Perform convolution to generate the projected image

[0152] Pixel weight calculation module 78 is used to use Calculate the projection image With the high-resolution particle image grayscale The ratio of each pixel

[0153] Update the spatial voxel grayscale module 79 to use the formula Update the grayscale of the spatial voxel. The grayscale of the voxel of the spatial particle field is updated with E(X j , Y j , Z j ) means, (X j , Y j , Z j ) is the coordinate of voxel j.

[0154] The iteration module 80 is used to repeatedly execute the particle projection matrix acquisition module, the projection image generation module, the pixel proportion calculation module and the space voxel grayscale update module until the preset number of iterations of the particle image is reached.

[0155] In some embodiments, the approximate weight matrix calculation module 75 is specifically configured to:

[0156] Assuming that each spatial voxel is projected to the pixel center, the approximate weight matrix is ​​calculated once:

[0157] The voxel simplified into a sphere is projected onto the high-resolution pixel plane, and the circle is called the voxel projection circle;

[0158] Calculate the area A of the high-resolution pixels overlapping the projected circle mn , m, n represent the indices of the overlapping pixels. If the pixels are completely within the projection circle, the overlapping area is d pixel 2 / R 2 If there is partial overlap, it can be approximated by the Monte Carlo method, where d pixel is the original pixel size, R represents the division level, and the optimal value of R is 4;

[0159] Calculate the approximate weight matrix The calculation formula is A mn / A projection , where A projection The area of ​​the voxel projection circle is represented by the formula Calculated;

[0160] According to the mapping function, the sphere centers of all non-zero voxels are projected onto the high-resolution pixel plane, the pixel centers closest to them are found, and the initial grayscale values ​​of these non-zero voxels are assigned to the corresponding pixels to obtain the particle projection matrix. If there are multiple spatial voxels corresponding to one pixel, the grayscale values ​​of these spatial voxels are accumulated and assigned to the pixel.

[0161] In some embodiments, the spatial voxel grayscale update module 79 is specifically configured to:

[0162] Residual Find the logarithm and get

[0163] Using the approximate weight matrix right Perform convolution operation to obtain

[0164] Exponentiate the result of the convolution operation to obtain Where r represents the pixel that overlaps with the projected circle of voxel j;

[0165] Will and Multiply and update the grayscale of the spatial voxel.

[0166] The convolution reconstruction system for tomographic particle image velocimetry of this embodiment has the same implementation process, method, and effect as the convolution reconstruction method for tomographic particle image velocimetry described in the first embodiment, and will not be described in detail here.

[0167] Example 3

[0168] like Figure 10 As shown, this embodiment relates to an electronic device including at least one processor and a memory communicatively connected to the at least one processor. The memory stores a computer program that can be run by the at least one processor, and the computer program is executed by the at least one processor so that the at least one processor can execute the convolution reconstruction method of a tomographic particle image velocimetry of Example 1, and realize the corresponding beneficial effects of the convolution reconstruction method of a tomographic particle image velocimetry, which will not be repeated here. The electronic device provided in this embodiment can be a personal computer, such as a desktop computer, an all-in-one computer, a laptop computer, a tablet computer, etc., and can also be a terminal device such as a mobile phone, a wearable device, a handheld computer, etc. The electronic device is merely an example and should not bring any limitations to the functions and scope of use of the embodiments of the present invention.

[0169] The components of the electronic device 3 may include, but are not limited to: the at least one processor 4 mentioned above, the at least one memory 5 mentioned above, and a bus 6 connecting different system components (including the memory 5 and the processor 4).

[0170] The bus 6 includes a data bus, an address bus, and a control bus.

[0171] The memory 5 may include a volatile memory, such as a random access memory (RAM) 51 and / or a cache memory 52 , and may further include a read-only memory (ROM) 53 .

[0172] The memory 5 may also include a program / utility 55 having a set (at least one) of program modules 54, such program modules 54 including but not limited to: an operating system, one or more application programs, other program modules and program data, each of which or some combination may include an implementation of a network environment.

[0173] The processor 4 executes various functional applications and data processing by running the computer program stored in the memory 5, such as the above-mentioned convolution reconstruction method for tomographic particle image velocimetry.

[0174] The electronic device 3 can also communicate with one or more external devices 7 (e.g., keyboard, pointing device, etc.). Such communication can be performed through an input / output (I / O) interface 8. Furthermore, the electronic device 3 can also communicate with one or more networks (e.g., a local area network (LAN), a wide area network (WAN), and / or a public network, such as the Internet) through a network adapter 9. Figure 5 As shown, the network adapter 9 communicates with other modules of the electronic device 3 via the bus 6. It should be understood that although Figure 5 Not shown, other hardware and / or software modules may be used in conjunction with the electronic device 3, including but not limited to: microcode, device drivers, redundant processors, external disk drive arrays, RAID (disk array) systems, tape drives, and data backup storage systems.

[0175] It should be noted that although several units / modules or sub-units / modules of the electronic device are mentioned in the detailed description above, this division is merely exemplary and not mandatory. In fact, according to embodiments of the present invention, the features and functions of two or more units / modules described above may be embodied in a single unit / module. Conversely, the features and functions of a single unit / module described above may be further divided and embodied by multiple units / modules.

[0176] Example 4

[0177] The present invention relates to a computer-readable storage medium, which stores instructions. When the instructions are executed, a convolution reconstruction method for tomographic particle image velocimetry is executed. The execution process, method and effect of the convolution reconstruction method for tomographic particle image velocimetry are the same as those described in Example 1, and will not be repeated here.

[0178] It should be noted that, in this document, the terms "comprises," "includes," or any other variations thereof are intended to encompass non-exclusive inclusion, such that a process, method, article, or apparatus comprising a series of elements includes not only those elements but also other elements not explicitly listed, or elements inherent to such process, method, article, or apparatus. In the absence of further limitations, an element defined by the phrase "comprising a ..." does not exclude the presence of other identical elements in the process, method, article, or apparatus comprising the element.

[0179] The above are only preferred embodiments of the present invention and are not intended to limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made using the contents of the present invention description and drawings, or directly or indirectly applied in other related technical fields, are also included in the patent protection scope of the present invention.

Claims

1. A convolution reconstruction method for tomographic particle image velocimetry, characterized in that: Used for reconstruction of three-dimensional flow fields, including: S101 obtains particle images required for reconstruction and sets reconstruction parameters, wherein the particle images are multiple particle images obtained by taking multiple cameras set at different angles; S102 initializes the discretized spatial voxels; S103: re-dividing the original pixel grid of the particle image based on the division level to obtain a high-resolution pixel grid, wherein the original pixel size in the high-resolution pixel grid is converted into an interpolated pixel size; S104 interpolates the particle image on the original pixel plane to the high-resolution pixel grid to obtain a high-resolution particle image grayscale image. S105 Calculate the approximate weight matrix of the spatial voxels according to the division level Specifically, it includes: assuming that each spatial voxel is projected to the pixel center, then calculating an approximate weight matrix, specifically: projecting the voxel simplified into a sphere onto the high-resolution pixel plane, the circle is called the voxel projection circle; calculating the area A of the high-resolution pixel overlapping with the projection circle mn , m, n represent the indices of the overlapping pixels. If the pixels are completely within the projection circle, the overlapping area is d pixel 2 / R 2 If there is partial overlap, it can be approximated by the Monte Carlo method, where d pixel is the original pixel size, R represents the division level, and the optimal value of R is 4; calculate the approximate weight matrix The calculation formula is A mn / A projection , where A projection The area of ​​the voxel projection circle is represented by the formula Calculation: According to the mapping function, the sphere centers of all non-zero voxels are projected onto the high-resolution pixel plane, the pixel centers closest to them are found, and the initial grayscale values ​​of these non-zero voxels are assigned to the corresponding pixels to obtain the particle projection matrix. If there are multiple voxels corresponding to one pixel, the grayscale values ​​of these voxels are accumulated and assigned to the pixel; S106 projects the sphere centers of all non-zero voxels onto the high-resolution pixel plane according to the mapping function, finds the pixel centers closest to them, assigns the initial grayscale values ​​of the non-zero voxels to the corresponding pixels, and obtains a particle projection matrix; S107 uses the approximate weight matrix as a convolution kernel and performs convolution with the particle projection matrix to generate a projection image S108 Exploitation Calculate the projection image With the high-resolution image The ratio of each pixel S109 Utilization Formula Update the grayscale of the spatial voxel, and the grayscale of the spatial voxel is expressed by E(X j , Y j , Z j ) means, (X j , Y j , Z j ) is the coordinate of the voxel in the j-space volume.

2. The convolution reconstruction method for tomographic particle image velocimetry according to claim 1, characterized in that: The multiple cameras are four cameras arranged in an X-shaped orthogonal distribution, and the angle between two adjacent cameras is 60 degrees. The reconstruction parameters include at least a mapping function coefficient, an iteration number, and a threshold coefficient. The method further includes: S110 repeatedly executing all steps S106 to S109 until a preset number of iterations of the particle image is reached.

3. The convolution reconstruction method for tomographic particle image velocimetry according to claim 1, characterized in that: The utilization formula Update the grayscale of the spatial voxel. The steps are: Residual Find the logarithm and get Using the approximate weight matrix right Perform convolution operation to obtain Exponentiate the result of the convolution operation to obtain Where r represents the pixel that overlaps with the projected circle of voxel j; Will and Multiply and update the grayscale of the spatial voxel.

4. The convolution reconstruction method for tomographic particle image velocimetry according to claim 1, characterized in that: The mapping function is: F(x)=a0+a1X+a2Y+a3X 2 +a4XY+ <h2 style=";text-align:left;direction:ltr">a5Y<h2 style=";text-align:left;direction:ltr"> 2 <h2 style=";text-align:left;direction:ltr"> +a6X<h2 style=";text-align:left;direction:ltr"> 3 <h2 style=";text-align:left;direction:ltr"> +a7X<h2 style=";text-align:left;direction:ltr"> 2 <h2 style=";text-align:left;direction:ltr"> Y+a8XY<h2 style=";text-align:left;direction:ltr"> 2 <h2 style=";text-align:left;direction:ltr"> +a9Y<h2 style=";text-align:left;direction:ltr"> 3 <h2 style=";text-align:left;direction:ltr">F(y)=b0+b1X+b2Y+b3X<h2 style=";text-align:left;direction:ltr"> 2 <h2 style=";text-align:left;direction:ltr"> +b4XY+ <h2 style=";text-align:left;direction:ltr">b5Y<h2 style=";text-align:left;direction:ltr"> 2 <h2 style=";text-align:left;direction:ltr"> +b6X<h2 style=";text-align:left;direction:ltr"> 3 <h2 style=";text-align:left;direction:ltr"> +b7X<h2 style=";text-align:left;direction:ltr"> 2 <h2 style=";text-align:left;direction:ltr"> Y+b8XY<h2 style=";text-align:left;direction:ltr"> 2 <h2 style=";text-align:left;direction:ltr"> +b9Y<h2 style=";text-align:left;direction:ltr"> 3 Where F(x) and F(y) represent the two-dimensional pixel coordinates, (X, Y) represents the three-dimensional space coordinates of the voxel, and a i and b i is the coefficient obtained by the least squares method, which is a function of the calibration position Z.

5. A convolution reconstruction system for tomographic particle image velocimetry, characterized in that: Used for reconstruction of three-dimensional flow fields, including: A reconstruction setting module is used to obtain particle images required for reconstruction and set reconstruction parameters, wherein the particle images are multiple particle images obtained by taking multiple cameras set at different angles; Discrete initialization module, used to initialize the discretized spatial voxels; a high-resolution pixel grid acquisition module, configured to re-divide the original pixel grid of the particle image based on the division levels to obtain a high-resolution pixel grid, wherein the original pixel size in the high-resolution pixel grid is converted into an interpolated pixel size; The interpolation module is used to interpolate the particle image on the original pixel plane to the high-resolution pixel grid to obtain the high-resolution particle image grayscale An approximate weight matrix calculation module is used to calculate the approximate weight matrix of the spatial voxels according to the division level The particle projection matrix acquisition module is used to project the sphere centers of all non-zero voxels onto the high-resolution pixel plane according to the mapping function, find the pixel centers closest to them, assign the initial grayscale values ​​of the non-zero voxels to the corresponding pixels, and obtain the particle projection matrix. A projection image generation module is used to use the approximate weight matrix As the convolution kernel and particle projection matrix Perform convolution to generate the projected image Pixel weight calculation module, used to use Calculate the projection image With the high-resolution particle image grayscale The ratio of each pixel Update the spatial voxel grayscale module to use the formula Update the grayscale of the spatial voxel. The grayscale of the voxel of the spatial particle field is updated with E(X j , Y j , Z j ) means, (X j ,Y j ,Z j ) is the coordinate of voxel j; The approximate weight matrix calculation module is specifically used to: Assuming that each spatial voxel is projected to the pixel center, the approximate weight matrix is ​​calculated once: the voxel simplified to a sphere is projected onto the high-resolution pixel plane. This circle is called the voxel projection circle. The area A of the high-resolution pixel overlapping with the projection circle is calculated. mn , m, n represent the indices of the overlapping pixels. If the pixels are completely within the projection circle, the overlapping area is d pixel 2 / R 2 If there is partial overlap, it can be approximated by the Monte Carlo method, where d pixel is the original pixel size, R represents the division level, and R is optimally 4; calculate the approximate weight matrix The calculation formula is A mn / A projection , where A projection The area of ​​the voxel projection circle is represented by the formula Calculation: According to the mapping function, the sphere centers of all non-zero voxels are projected onto the high-resolution pixel plane, the pixel centers closest to them are found, and the initial grayscale values ​​of these non-zero voxels are assigned to the corresponding pixels to obtain the particle projection matrix. If there are multiple spatial voxels corresponding to one pixel, the grayscale values ​​of these spatial voxels are accumulated and assigned to the pixel.

6. The convolution reconstruction system for tomographic particle image velocimetry according to claim 5, characterized in that: The multiple cameras are four cameras distributed orthogonally in an X shape, and the angle between two adjacent cameras is 60 degrees. The reconstruction parameters include at least a mapping function coefficient, an iteration number, and a threshold coefficient; It also includes: an iteration module for repeatedly executing the particle projection matrix acquisition module, the projection image generation module, the pixel proportion calculation module and the space voxel grayscale update module until the preset number of iterations of the particle image is reached.

7. The convolution reconstruction system for tomographic particle image velocimetry according to claim 5, characterized in that: The updating spatial voxel grayscale module is specifically used to: Residual Find the logarithm and get Using the approximate weight matrix right Perform convolution operation to obtain Exponentiate the result of the convolution operation to obtain Where r represents the pixel that overlaps with the projected circle of voxel j; Will and Multiply and update the grayscale of the spatial voxel.

8. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the computer program, the convolution reconstruction method for tomographic particle image velocimetry according to any one of claims 1 to 4 is implemented.

9. A computer-readable storage medium, characterized in that: The storage medium stores instructions, which, when executed, execute a convolution reconstruction method for tomographic particle image velocimetry according to any one of claims 1 to 4.

Images