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

Abstract

The invention relates to a convolution reconstruction method, system and device for tomography particle image velocimetry and a medium, and the method comprises the steps: S101, obtaining a particle image needed by reconstruction, and setting reconstruction parameters; s102, carrying out initialization on the discrete space voxels; s103, re-grading the original pixel grid of the particle image to obtain a high-resolution pixel grid; s104, obtaining the gray level of the high-resolution particle image; s105, calculating an approximate weight matrix of spatial voxels according to the division grades; s106, projecting sphere centers of all non-zero voxels to a high-resolution pixel plane according to a mapping function, and assigning a particle projection matrix by using a nearest pixel center; s107, using the approximate weight matrix as a convolution kernel to perform convolution with the particle projection matrix to generate a projection image; s108, calculating a ratio corresponding to each pixel between the projection image and the high-resolution image; and S109, updating the gray scale of the space voxel. According to the invention, the particle field reconstruction efficiency and quality of tomography particle image velocimetry are improved.

Description

Technical Field

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

[0002] Tomographic Particle Image Velocimetry (Tomo-PIV for short) technology is an important means for quantitatively measuring three-dimensional flow fields in experimental fluid mechanics. This technology records the particle motion in the flow field from multiple perspectives, uses advanced reconstruction algorithms to construct the three-dimensional particle distribution in the flow field, and finally calculates the three-dimensional velocity field distribution of the flow field through three-dimensional cross-correlation algorithms. The birth of Tomo-PIV technology not only marks the leap of flow measurement technology from two dimensions to three dimensions, but also provides a more powerful tool for studying complex flow phenomena such as turbulence and eddy currents.

[0003] The core step of Tomo-PIV technology is spatial gray reconstruction (also known as particle field reconstruction). The so-called spatial gray reconstruction is to find the spatial gray distribution E given the particle image P and the weight function W, which is the inverse process of the projection imaging WE = P. Since the number of discrete voxels is much larger than the number of known pixels, this inverse problem has an indefinite solution. The most common current reconstruction algorithm is the multiplicative algebraic reconstruction technique (MART for short), and its iterative formula is:

[0004]

[0005] where (x i , y i ) is the coordinate of the i-th pixel, P(x i , y i ) is the pixel gray value of the i-th pixel, (X j , Y j , Z j ) is the coordinate of the j-th voxel, B(X j , Y j , Z j ) is the voxel gray value of the j-th voxel, w i，j is the weight coefficient, indicating the contribution value of the j-th voxel to the pixel gray value of the i-th pixel. This contribution is related to the projection distance, and the smaller the distance, the greater the contribution. The range of the weight coefficient is between 0 and 1. k represents the number of iterations, and μ is the relaxation coefficient.

[0006] Although the MART method can effectively achieve three-dimensional reconstruction of the particle field, its greatest limitation lies in its low computational efficiency. The main reasons are as follows: First, in the calculation process of the MART method, the weight coefficient of each voxel is calculated separately, and the amount of calculation is proportional to the cube of the size of the reconstructed volume. Second, the MART method requires a large amount of memory and has high requirements for computing devices. These problems seriously affect the practical application of Tomo-PIV.

[0007] In view of the limitations of the MART method, scholars have conducted continuous research and proposed a variety of improved fast reconstruction algorithms. Among them, the multiplicative line of sight-simultaneous MART (MLOS-SMART) can significantly reduce memory occupancy while improving the reconstruction speed. Therefore, it is the most widely used fast reconstruction algorithm at present. However, this algorithm still does not increase the reconstruction speed to the fastest, and it 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 particle field reconstruction in Tomo-PIV is an urgent problem to be solved. Summary of the Invention

[0008] The present invention provides a convolutional 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 object, in the first aspect, the present invention provides a convolutional reconstruction method for tomographic particle image velocimetry, including:

[0010] S101 Obtain the particle images required for reconstruction and set the reconstruction parameters, where the particle images are multiple particle images taken by multiple cameras set at different angles;

[0011] S102 Initialize the discretized spatial voxels;

[0012] S103 Re-divide the original pixel grid of the particle image based on the division level to obtain a high-resolution pixel grid, and the original pixel size in the high-resolution pixel grid becomes the interpolated pixel size;

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

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

[0015] S106 Project the centers of all non-zero voxels onto the high-resolution pixel plane according to the mapping function, find the pixel centers closest to them, and assign the initial gray values of the non-zero voxels to the corresponding pixels to obtain a particle projection matrix;

[0016] S107 Use the approximate weight matrix as a convolution kernel to convolve with the particle projection matrix to generate a projection image

[0017] S108 Use to calculate the projection image and the ratio corresponding to each pixel between the projection image and the high-resolution image

[0018] S109 Use the formula to update the gray value of the spatial voxel. The gray value of the spatial voxel is represented by E(X j , Y j , Z j ). E(X j , Y j , Z j ) is the coordinate of the voxel in the j space body.

[0019] Preferably, the multiple cameras are four cameras orthogonally distributed in an X shape, and the included angle between two adjacent cameras is 60 degrees. The reconstruction parameters at least include mapping function coefficients, the number of iterations, and threshold coefficients;

[0020] It further includes: S110 Repeat all steps of S106 - S109 until the preset number of iterations of the particle image is reached.

[0021] Preferably, the division level is 4.

[0022] Preferably, calculating the approximate weight matrix of the voxel according to the division level specifically includes: Assuming that each spatial voxel is projected onto the pixel center, calculate the approximate weight matrix once. Specifically:

[0023] Project the voxel simplified as a sphere onto the high-resolution pixel plane. This circle is called the voxel projection circle;

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

[0025] Calculate the approximate weight matrix The calculation formula is A mn / A projection , where A projectoin represents the area of the projection circle of the spatial voxel, which is calculated by the formula ;

[0026] Project the centers of the spheres of all non-zero voxels onto the high-resolution pixel plane according to the mapping function, find the pixel centers closest to them, and assign the initial gray values of these non-zero voxels to the corresponding pixels to obtain the particle projection matrix If multiple voxels correspond to one pixel, accumulate the gray values of these voxels and assign them to the pixel.

[0027] Preferably, using the formula to update the gray value of the spatial voxel, the steps are as follows:

[0028] Take the logarithm of the residual to obtain

[0029] Perform a convolution operation on using the approximate weight matrix to obtain

[0030] Take the exponential of the result obtained from the convolution operation to obtain where r represents the pixel overlapping with the projection circle of the j-th voxel;

[0031] Multiply by to update the gray value of the spatial voxel.

[0032] Preferably, the mapping function is:

[0033] F(x) = a 0 + a 1 X + a 2 Y + a 3 X 2 + a 4 XY + a 5 Y 2 + a 6 X 3 + a 7 X 2 Y + a 8 XY 2 + a 9 Y 3

[0034]

[0035] To achieve the above object, in a second aspect, for the reconstruction of a three-dimensional flow field, it includes:

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

[0037] A discrete initialization module, configured 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 a division level to obtain a high-resolution pixel grid, and the original pixel size in the high-resolution pixel grid becomes an interpolated pixel size;

[0039] An interpolation module, configured to interpolate the particle image on the original pixel plane onto the high-resolution pixel grid to obtain the gray level of the high-resolution particle image

[0040] An approximate weight matrix calculation module, configured to calculate an approximate weight matrix of the spatial voxels according to the division level

[0041] A particle projection matrix acquisition module, configured to project the centers of all non-zero voxels onto the high-resolution pixel plane according to the mapping function, find the pixel centers closest to them, and assign the initial gray level values of the non-zero voxels to the corresponding pixels to obtain a particle projection matrix

[0042] A projection image generation module, configured to use the approximate weight matrix as a convolution kernel and the particle projection matrix to perform convolution to generate a projection image

[0043] A pixel ratio calculation module, configured to use to calculate the ratio corresponding to each pixel between the projection image and the gray level of the high-resolution particle image

[0044] An updated spatial voxel gray level module, configured to update the gray level of the spatial voxels using the formula to update the gray level of the spatial voxels, and the gray level of the voxels in the spatial particle field is represented by E(X j , Y j , Z j ), and (X j , Y j , Z j ) is the coordinate of the j-th voxel.

[0045] ​Preferably, the multiple cameras are four cameras orthogonally distributed in an X shape, and the included angle between two adjacent cameras is 60 degrees. The reconstruction parameters at least include mapping function coefficients, the number of iterations, and threshold coefficients. It further includes: an iteration module, configured to repeatedly execute the particle projection matrix obtaining module, the projection image generating module, the pixel ratio calculating module, and the spatial voxel gray level updating module until the preset number of iterations of the particle image is reached.

[0046] Preferably, the approximate weight matrix calculating module is specifically configured to:

[0047] Assume that each spatial voxel projects onto the pixel center, and calculate an approximate weight matrix once:

[0048] Project the voxel simplified as a sphere onto a high-resolution pixel plane, and this circle is called the voxel projection circle;

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

[0050] Calculate the approximate weight matrix The calculation formula is A mn / A projection where A projection represents the area of the spatial voxel projection circle, which is calculated by the formula ;

[0051] According to the mapping function, project the centers of all non-zero voxels onto the high-resolution pixel plane, find the pixel centers closest to them, and assign the initial gray levels of these non-zero voxels to the corresponding pixels to obtain the particle projection matrix If multiple spatial voxels correspond to one pixel, the gray levels of these spatial voxels are accumulated and assigned to this pixel.

[0052] Preferably, the spatial voxel gray level updating module is specifically configured to:

[0053] Take the logarithm of the residual to obtain

[0054] Convolve with the approximate weight matrix to obtain

[0055] Take the exponent of the result obtained from the convolution operation to obtain Among them, r represents the pixels overlapping with the projection circle of the j - voxel;

[0056] Multiply with to update the gray value of the spatial voxel.

[0057] To achieve the above object, in a third aspect, the present invention further relates to an electronic device, including a memory, a processor, and a computer program stored on the memory and executable on the processor, wherein the processor implements the above - mentioned tomographic particle image velocimetry convolution reconstruction method when executing the computer program.

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

[0059] A tomographic particle image velocimetry convolution reconstruction method, system, device, and medium according to the present invention have the following beneficial effects compared with the prior art:

[0060] Based on the MLOS - SMART method, a Convolutional MART (Conv - MART for short) algorithm is proposed, aiming to improve the particle field reconstruction efficiency of Tomo - PIV. This algorithm re - partitions the pixel grid to obtain high - resolution pixels, uses an approximate weight matrix to replace different weight matrices to reduce the calculation time, and introduces a convolution operation on the basis of the SMART method to further accelerate the calculation process, resulting in the Conv - MART algorithm. Compared with the MLOS - SMART method, the advantages of this algorithm are: 1) It reduces the calculation time of the weight matrix and has fewer calculation steps, greatly improving the reconstruction efficiency; 2) When the particle image contains noise or the particle concentration is high, the Conv - MART method has higher reconstruction quality than the MLOS - SMART method. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 is a flowchart of a tomographic particle image velocimetry convolution reconstruction method in Embodiment 1 of the present invention;

[0062] Figure 2 is a schematic diagram of dividing a pixel grid of a tomographic particle image velocimetry convolution reconstruction method in Embodiment 1 of the present invention;

[0063] Figure 3 is a flowchart of step 105 of a tomographic particle image velocimetry convolution reconstruction method in Embodiment 1 of the present invention;

[0064] Figure 4It is the flowchart of step 109 of a convolution reconstruction method for tomographic particle image velocimetry in Embodiment 1 of the present invention;

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

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

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

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

[0069] Figure 9 It is the structural schematic diagram of a convolution reconstruction system for tomographic particle image velocimetry in Embodiment 2 of the present invention;

[0070] Figure 10 It is the structural schematic diagram of an electronic device in Embodiment 3 of the present invention. Specific embodiments

[0071] The present invention will be further described in detail below with reference to the drawings and embodiments. It can be understood that the specific embodiments described herein are only used to explain the present invention, rather than limiting the present invention. Additionally, it should be noted that for the convenience of description, only parts related to the present invention are shown in the drawings, rather than all the structures.

[0072] Embodiment 1

[0073] A convolution reconstruction method for tomographic particle image velocimetry. Please refer to Figures 1-8 , which is used for the reconstruction of three-dimensional flow fields. The present invention proposes a convolution multiplicative algebraic reconstruction method (Conv-MART) aiming to improve the particle field reconstruction efficiency of Tomo-PIV. In the present invention, the pixel gray level of the particle image is represented by P(x i , y i ), where (x i , y i ) are the coordinates of the pixel. The voxel gray level of the spatial particle field is represented by E(X j , Y j , Z j ), where (X j , Y j , Z j ) are the coordinates of the j voxel. The weight coefficient between the pixel and the voxel is w i，jThe voxels are simplified to spheres, and the pixels remain square. The present invention mainly re - divides the particle image to obtain 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 calculation to obtain the reconstructed particle field. In addition, the present invention can conveniently use a graphics processing unit (GPU) to further accelerate the reconstruction.

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

[0075] S101: Obtain the particle image required for reconstruction and set the reconstruction parameters, where the particle image is a plurality of particle images taken by multiple cameras set at different angles.

[0076] In this embodiment, the multiple cameras are four cameras orthogonally distributed in an X - shape, and the included angle between two adjacent cameras is 60 degrees. The reconstruction parameters at least include the mapping function coefficient, the number of iterations, and the threshold coefficient. If a GPU needs to be combined, the particle image and the mapping function coefficient in this step need to be initialized as GPU data.

[0077] S102: Initialize the discretized spatial voxels.

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

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

[0080] Among them, the pixel size changes from the original d pixel (generally 1 pixel) to d pixel / R, where R represents the division level (Resize Level). Figure 2 FIG. is a schematic diagram of dividing the pixel grid. The divided pixel size is called the interpolated pixel size. The present invention has been tested and found that the value of R is most suitable when taken as 4, and at this time the interpolated pixel size is d pixel / 4.

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

[0082] Among them, the interpolation method includes but is not limited to bilinear interpolation, polynomial interpolation, nearest - neighbor interpolation, etc.

[0083] S105: Calculate the approximate weight matrix of the spatial voxels according to the divided levels

[0084] In this step, it is assumed that each spatial voxel is projected onto the pixel center. Without the projection process, the approximate weight matrices calculated for different voxels are the same. Therefore, it only needs to be calculated once. The higher the resolution of the generated image, the more valid this assumption is. Figure 3 This is the flowchart of this step. The specific process includes S1051 - S1053:

[0085] Assume that each spatial voxel is projected onto the pixel center, then calculate the approximate weight matrix once. Specifically:

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

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

[0088] S1053: Calculate the approximate weight matrix The calculation formula is A mn / A projection , where A projection represents the area of the projection circle of the spatial voxel, which is calculated by the formula

[0089] The matrix is a square matrix. The element values and matrix dimensions are related to the ratio M of the divided level R and d voxel and d pixel .

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

[0091] In this embodiment, the mapping function used is formula (1, 2). The pinhole imaging model or other relationships can also be used.

[0092] F(x) = a​0 +a 1 X + a 2 Y + a 3 X 2 + a 4 XY + a 5 Y 2 + a 6 X 3 + a 7 X 2 Y + a 8 XY 2 + a 9 Y 3 (1)

[0093] F(y) = b 0 + b 1 X + b 2 Y + b 3 X 2 + b 4 XY + b 5 Y 2 + b 6 X 3 + b 7 X 2 Y + b 8 XY 2 + b 9 Y 3 (2).

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

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

[0096]

[0097] Among them, "*" represents the convolution operation on the matrix The convolution kernel is the weight matrix obtained in step 105 Use the convolution formula to generate a projection image The calculation formula for this step is as follows:

[0098]

[0099] Among them, "*" represents the convolution operation on the matrix The convolution kernel is the weight matrix obtained in S105

[0100] S108: Use to calculate the projection image The ratio corresponding to each pixel between and the high-resolution image

[0101] S109: Update the gray value of the spatial voxel using formula (3), where the gray value of the spatial voxel is represented by E(X j , Y j , Z j ), and (X j , Y j , Z j ) is the coordinate of the voxel in the j-th space.

[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: Take the logarithm of the residual to obtain

[0106] S1092: Perform a convolution operation on using the approximate weight matrix to obtain

[0107] S1093: Take the exponent of the result obtained from the convolution operation to obtain where r represents the pixels overlapping the projection circle of the j-th voxel.

[0108] It can be mathematically proven that is equivalent to the of the SMART method and has higher computational efficiency. Here, r represents the pixels overlapping the projection circle of the j-th voxel.

[0109] S1094: Multiply by to update the gray value of the spatial voxel.

[0110] S110: Repeat all the above steps of 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 but only simple algebraic operations, and can be fully implemented on the GPU, thereby further improving the reconstruction efficiency. Through the above steps, that is, on the premise of ensuring the reconstruction quality, the reconstructed particle field can be calculated in a very short time.

[0112] 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 simulation. In this example, a vortex ring is simulated, with the number of voxels of the vortex ring being 700×700×140 voxels and the resolution being 0.05 mm / voxel. 24,500 tracer particles are randomly distributed therein, and the gray scale of each spatial particle follows a Gaussian distribution. Then, the particles in the vortex ring are projected onto four cameras to obtain four particle images. These four cameras are distributed in an X shape, and the included angle between two adjacent cameras is 60°. The resolution of the particle image is 900×900 pixels, and the size of the pixel is 0.05 mm. We use the Conv-MART method and the MLOS-SMART method respectively to reconstruct the particle fields at two consecutive moments in the vortex ring, and test 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.05 ppp (particle per pixel) in the range of 700×700×140 voxels, which is the particle field at time A. Then, calculate the displacement of the particles after a short time according to 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 included angle between the lines of sight of adjacent cameras being 60 degrees.

[0117] Step 3: Use the MLOS method to initialize the gray scale 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 / 4 mm, re-divide the pixel plane, and obtain a high-resolution pixel plane.

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

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

[0121] Step 7: Use the mapping functions (1) and (2) to calculate and obtain the particle projection matrix Then use the convolution formula to generate a projection image

[0122] Step 8: Use the formula to calculate the ratio corresponding to each pixel between the projected image and the high-resolution image .

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

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

[0125] Step 11: Adopt 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 comparison of the final calculation time results is as Figure 5 shown, which respectively shows the average time for the MLOS-SMART method, the Conv-MART method, and the GPU-Conv-MART method to reconstruct a particle field.

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

[0127] There are only differences in the reconstruction time between the Conv-MART method and the GPU-Conv-MART method, and there are no differences in the reconstructed particle field and the calculated velocity field. Therefore, Figure 6 only the velocity field of the Conv-MART method is shown. The arrows in the figure represent the solved velocity vectors, and the purple isosurface is the Q-criterion vorticity isosurface.

[0128] Example 2:

[0129] Example 2 is a left ventricular Tomo - PIV experiment, mainly used to measure the blood flow structure in the left ventricle after replacing the native mitral valve with a mechanical valve. In this experiment, four high - speed cameras equipped with 100mm / f2.8 lenses are arranged in a straight line, with a sampling frequency of 200Hz, and the angle between two adjacent cameras is 36°. The resolution of the cameras is 2016×2016, and the pixel size is 0.05mm. The tracer particles used in the experiment are rhodamine B fluorescent particles with a diameter of 50μm, and the particle concentration is approximately 0.045 ppp. In addition, in order to simulate the blood flow state in the left ventricle, the solution used in the experiment is a water - glycerol mixture, with 62% water and 38% glycerol by volume. The refractive index of this solution is 1.395, and its viscosity is 3.86 mPa·s at 25.6°C. The volume of the reconstructed measurement volume is taken as 65×70×57mm3, and the number of voxels included is 1300×1400×1100. In this example, we use these four high - speed cameras to capture a series of consecutive particle images in time series, select the particle images at the first two moments (set as moments A and B), and use the MLOS - SMART method and the Conv - MART method for reconstruction.

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

[0131] Step 2: Initialize the gray level of the particle field using the MLOS method, and find all non - zero voxels.

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

[0133] Step 4: Interpolate the particle image onto the high - resolution pixel plane using the bilinear interpolation method to obtain

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

[0135] Step 6: Calculate the particle projection matrix using the mapping functions (1) and (2) Then use the convolution formula to generate the projection image

[0136] Step 7: Use the formula to calculate the projection image and the high - resolution image The ratio corresponding to each pixel therebetween.

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

[0138] Step Nine: Input the particle image at time B, keep the other reconstruction parameters unchanged, and repeat Steps One to Eight to reconstruct the particle field at time B.

[0139] Step Ten: Adopt the GPU version of the algorithm, repeat Steps One to Nine to reconstruct the particle fields at times A and B. This step is used to test the acceleration performance of the Conv-MART method combined with GPU. The comparison results of the final calculation time are as Figure 7 shown, which respectively show the average time for the MLOS-SMART method, the Conv-MART method, and the GPU-Conv-MART method to reconstruct a particle field.

[0140] Step Eleven: Calculate the particle field reconstructed by the Conv-MART method using the multi-layer iterative volume deformation method to obtain the velocity field, and the result is as Figure 8 shown.

[0141] There are only differences in the reconstruction time between the Conv-MART method and the GPU-Conv-MART method, and there are no differences in the reconstructed particle field and the calculated velocity field. Therefore, Figure 8 only the velocity field of the Conv-MART method is shown. The arrows in the figure represent the solved velocity vectors, and the green isosurface is the Q-criterion vorticity isosurface.

[0142] S10: Example Two

[0143] A convolutional reconstruction system for tomographic particle image velocimetry, for reference please see Figures 1-8, for the reconstruction of three-dimensional flow fields, the present invention proposes a Convolution Multiplication Algebraic Reconstruction Technique (Conv-MART), aiming to improve the particle field reconstruction efficiency of Tomo-PIV. In this embodiment, the voxels are simplified to spheres, and the pixels are still squares. This system mainly re-divides the particle images to obtain 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 the reconstructed particle field. In addition, this system can easily utilize a Graphics 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, a server, an intelligent terminal, a 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 updated spatial voxel grayscale module 79, and an iterative module 80.

[0144] The reconstruction setting module 71 is used to obtain the particle images required for reconstruction and set the reconstruction parameters. Among them, the particle images are multiple particle images captured by multiple cameras with different angular settings.

[0145] The multiple cameras are four cameras arranged in an X-shaped orthogonal distribution, and the included angle between two adjacent cameras is 60 degrees. The reconstruction parameters at least include the mapping function coefficient, the number of iterations, and the 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 images based on the division level to obtain a high-resolution pixel grid, and the original pixel size in the high-resolution pixel grid becomes the interpolated pixel size.

[0148] The interpolation module 74 is used to interpolate the particle images on the original pixel plane onto 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 centers of all non-zero voxels onto the high-resolution pixel plane according to the mapping function, find the nearest pixel centers to them, and assign the initial grayscale values of the non-zero voxels to the corresponding pixels to obtain the particle projection matrix

[0151] A projection image generation module 77, configured to use the approximate weight matrix as a convolution kernel and a particle projection matrix to perform convolution to generate a projection image

[0152] A pixel ratio calculation module 78, configured to use to calculate the ratio corresponding to each pixel between the projection image and the gray level of the high-resolution particle image

[0153] An updated spatial voxel gray level module 79, configured to use a formula to update the gray level of the spatial voxel. The gray level of the voxel of the spatial particle field is represented by E(X j j , Y j , Z j ), (X j , Y j , Z j ) being the coordinates of the j voxel.

[0154] An iteration module 80, configured to repeatedly execute the particle projection matrix acquisition module, the projection image generation module, the pixel ratio calculation module, and the updated spatial voxel gray level 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] Assume that each spatial voxel projects onto the pixel center, and calculate the approximate weight matrix once:

[0157] Project the voxel simplified as a sphere onto the high-resolution pixel plane, and this circle is called the voxel projection circle;

[0158] Calculate the area A of the high-resolution pixels overlapping with the projection circle mn , where m and n represent the indices of the overlapping pixels. If the pixel is completely within the projection circle, the overlapping area is d pixel 2 / R 2 ; if there is partial overlap, it can be approximately calculated by the Monte Carlo method, where d pixel is the original pixel size, and 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 represents the area of the spatial voxel projection circle, which is calculated by the formula ;

[0160] Project the centers of all non-zero voxels onto the high-resolution pixel plane according to the mapping function, find the pixel centers closest to them, and assign the initial gray values of these non-zero voxels to the corresponding pixels to obtain the particle projection matrix. If multiple spatial voxels correspond to one pixel, the gray values of these spatial voxels are accumulated and assigned to the pixel.

[0161] In some embodiments, the spatial voxel gray module 79 is updated, specifically for:

[0162] Take the logarithm of the residual to obtain

[0163] Use the approximate weight matrix to perform a convolution operation on to obtain

[0164] Take the exponent of the result obtained from the convolution operation to obtain where r represents the pixels overlapping with the projection circle of the j-th voxel;

[0165] Multiply by to update the gray value of the spatial voxel.

[0166] The implementation process, method, and effect of the tomographic particle image velocimetry convolution reconstruction system in this embodiment are the same as those of the tomographic particle image velocimetry convolution reconstruction method described in Embodiment 1, and will not be elaborated here.

[0167] Embodiment 3

[0168] As Figure 10 shown, this embodiment relates to an electronic device including at least one processor and a memory communicatively connected to the at least one processor. Among them, 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 tomographic particle image velocimetry convolution reconstruction method of Embodiment 1 and achieve the corresponding beneficial effects of the tomographic particle image velocimetry convolution reconstruction method, which will not be elaborated 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 only an example and should not bring any limitations to the functions and usage ranges of the embodiments of the present invention.

[0169] The components of the electronic device 3 may include, but are not limited to: at least one of the above-mentioned processors 4, at least one of the above-mentioned memories 5, 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 volatile memory, such as random access memory (RAM) 51 and / or cache memory 52, and may further include read-only memory (ROM) 53.

[0172] The memory 5 may further include a program / utilities 55 having a set (at least one) of program modules 54. Such program modules 54 include, but are not limited to: an operating system, one or more application programs, other program modules, and program data. Each or some combination of these examples may include the implementation of a network environment.

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

[0174] The electronic device 3 may also communicate with one or more external devices 7 (such as a keyboard, a pointing device, etc.). Such communication may be carried out through an input / output (I / O) interface 8. And, the electronic device 3 may also communicate with one or more networks (such as a local area network (LAN), a wide area network (WAN), and / or a public network, such as the Internet) through a network adapter 9. As Figure 5 shown, the network adapter 9 communicates with other modules of the electronic device 3 through the bus 6. It should be understood that although Figure 5 not shown in the figure, other hardware and / or software modules may be used in combination 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, etc.

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

[0176] Embodiment 4

[0177] The present invention relates to a computer-readable storage medium storing instructions which, when running, execute a convolution reconstruction method for tomographic particle image velocimetry in Embodiment 1. The implementation process, method, and effects during its operation are the same as those of the convolution reconstruction method for tomographic particle image velocimetry described in Embodiment 1, and will not be elaborated herein.

[0178] It should be noted that, in this text, the term "including", "comprising", or any other variant thereof is intended to cover non-exclusive inclusion, such that a process, method, article, or device including a series of elements not only includes those elements but also includes other elements not explicitly listed, or further includes elements inherent to such process, method, article, or device. Without further limitation, an element defined by the statement "including one..." does not exclude the existence of another identical element in the process, method, article, or device including that element.

[0179] The above are only the preferred embodiments of the present invention and do not limit the patent scope of the present invention accordingly. Any equivalent structural or equivalent process transformation made by using the specification and drawings of the present invention, or directly or indirectly applied in other related technical fields, shall be equally 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 pictures with 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 changed 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 S106: projecting the sphere centers of all non-zero voxels onto the high-resolution pixel plane according to the mapping function, finding the pixel centers closest to them, assigning the initial grayscale values ​​of the non-zero voxels to the corresponding pixels, and obtaining a particle projection matrix; S107: Use the approximate weight matrix as a convolution kernel and the particle projection matrix to perform convolution to generate a projection image. S108 Utilization 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 ) are the coordinates 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, the angle between two adjacent cameras is 60 degrees, and the reconstruction parameters at least include a mapping function coefficient, an iteration number, and a threshold coefficient; The method further comprises: S110 repeatedly executing all steps S106-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 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, which is: The voxel simplified to a sphere is projected onto the high-resolution pixel plane, and the circle is called the voxel projection circle; 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; Calculate the approximate weight matrix The calculation formula is A mn / A projection , where A projection The area of ​​the voxel projection circle in the space is represented by the formula Calculated; 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.

4. 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 we get Using the approximate weight matrix right Perform convolution operation to obtain Exponentially taking the result of the convolution operation, we get Where r represents the pixel that overlaps with the projection circle of voxel j; Will and Multiply and update the grayscale of the voxel in the space.

5. The convolution reconstruction method for tomographic particle image velocimetry according to claim 3, characterized in that: The mapping function is: <h2 style=";text-align:left;direction:ltr">F(x) = a0 + a1X + a2Y + a3X<h2 style=";text-align:left;direction:ltr"> 2 <h2 style=";text-align:left;direction:ltr"> +a4XY+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+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 spatial 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.

6. 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, 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; Discrete initialization module, used to initialize the discretized spatial voxels; A high-resolution pixel grid acquisition module, used for re-dividing 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 is changed into an interpolated pixel size; An interpolation module is used to interpolate the particle image on the original pixel plane to the high-resolution pixel grid to obtain a high-resolution particle image grayscale image. 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 spherical 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 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 in the spatial particle field is updated by E(X j , Y j , Z j ) means, (X j , Y j , Z j ) are the coordinates of voxel j.

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

8. The convolution reconstruction system for tomographic particle image velocimetry according to claim 1, characterized in that: The approximate weight matrix calculation module is specifically used for: 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, and the circle is called the voxel projection circle; 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; Calculate the approximate weight matrix The calculation formula is A mn / A projection , where A projection The area of ​​the voxel projection circle in the space is represented by the formula Calculated; 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.

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

10. 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 according to any one of claims 1 to 5 is implemented.

11. 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 5.

Images