Particle image segmentation method for particle image velocity measurement

Abstract

The invention relates to the technical field of image processing, solves the technical problem that particle information in an image is difficult to retain to the greatest extent while interference factors in shooting data are eliminated in a traditional method, and particularly relates to a particle image segmentation method for particle image velocity measurement. Comprising particle image acquisition, threshold segmentation, contour extraction and area statistics, non-boundary abnormal region filtering, background image generation and particle head portrait segmentation. According to the method, interference can be removed, and meanwhile particle information in the image can be reserved to the maximum extent. Interference factors in the particle image are identified and removed, and only effective particles are reserved; the method is high in universality, can provide more accurate and meticulous data support for subsequent research and analysis of the flow field, and has a wide application prospect.

Description

Technical Field

[0001] The present invention relates to the technical field of image processing, and in particular to a particle image segmentation method for particle image velocimetry. Background Art

[0002] In the technical field of image processing, particle image velocimetry (PIV) is an important technique for measuring fluid motion. Due to its characteristics of full flow field, high precision and no interference with fluid motion, it is widely used in the measurement of fluid velocity fields. Tracer particles are seeded in the area to be measured, and after illuminating the tracer particles with a sheet light source, a high-speed camera is used for analysis.

[0003] However, since particle image velocimetry uses image cross-correlation for velocity measurement, when there are interference factors (such as background objects, other light sources, etc.) in the photographed area, the measurement results will have large errors. Most existing methods eliminate the influence of background objects by setting a dynamic mask. Setting a mask means identifying the shape of the background objects in the image and generating a mask, and shielding the corresponding mask area in the particle image.

[0004] The existing methods mainly perform threshold segmentation directly after preprocessing the image with background objects, or train a segmentation model through deep learning. The former has a relatively single threshold setting, and the latter requires a large amount of materials for training, and the effect is not ideal for various complex and changeable background scenes. For example, a Chinese invention patent with the publication number CN115170639A and the invention name of an image processing method for three-dimensional particle image velocimetry of a robotic fish underwater. In addition, when setting a mask, the particles in this area will also be shielded, and the method of setting a mask cannot cope with the interference of background light sources. Therefore, accurately identifying the particle area and eliminating interference factors is of great significance for improving the accuracy and robustness of PIV measurement. Summary of the Invention

[0005] Aiming at the deficiencies of the prior art, the present invention provides a particle image segmentation method for particle image velocimetry, which solves the technical problem that it is difficult for traditional methods to eliminate interference factors in the photographed data while maximizing the retention of particle information in the image.

[0006] To solve the above technical problems, the present invention provides the following technical solutions: A particle image segmentation method for particle image velocimetry, the method includes the following steps:

[0007] S1. Obtain a particle image I located in the measurement area through a two-dimensional PIV measurement device;

[0008] S2. Define a global segmentation threshold T and a local segmentation threshold T *, and segment the particle image I to obtain a binary image bw1 with the foreground and background separated, and a binary image bw2 that segments the particles overlapping with the background in the particle image I;

[0009] S3. Perform connected component calculation on the binary image bw1 and the binary image bw2, and remove extreme values to determine the particle reference area S O ;

[0010] S4. Use the particle reference area S O to filter the connected components that are non-boundary abnormal regions to obtain a new binary image bw;

[0011] S5. Generate a background image B that is consistent with the grayscale value of the particle image I based on the binary image bw;

[0012] S6. Segment the particle image I according to the background image B to obtain the particle image I * .

[0013] Further, the process of obtaining the particle image I includes:

[0014] Build a particle image velocimetry environment using a high-speed camera, a plane laser, and a water tank, keep the water flow in the water tank in a flowing state, then adjust the plane laser settings to obtain a measurement plane, subsequently sprinkle tracer particles to form a flow field environment to be measured, and finally collect the particle image I in the measurement area through the high-speed camera.

[0015] Further, in step S2, the specific process includes the following steps:

[0016] S21. Determine the optimal threshold T i as the global segmentation threshold T;

[0017] S22. Perform binarization processing on the particle image I using the global segmentation threshold T to obtain a binary image bw1. The expression of the pixel bw1(x, y) in the x-th column and y-th row of the binary image bw1 is:

[0018]

[0019] S23. Determine the optimal threshold T i * as the local segmentation threshold T * ;

[0020] S24. Perform binarization processing on the particle image I using the local segmentation threshold T * to obtain a binary image bw2. The expression of the pixel bw2(x, y) in the x-th column and y-th row of the binary image bw2 is:

[0021]

[0022] In the formula, I(x, y) represents any pixel point in the particle image I.

[0023] Furthermore, in step S21, the specific process includes the following steps:

[0024] S211. Calculate the probability P(i) of each gray level i according to the gray level histogram of the particle image I. The calculation formula is:

[0025]

[0026] In the formula, n i is the number of pixels of gray level i; N is the total number of pixels in the particle image I;

[0027] S212. Calculate the foreground mean μ 0 (T) and the background mean μ 1 (T). The calculation formulas for the foreground mean μ 0 (T) and the background mean μ 1 (T) are respectively:

[0028]

[0029] In the formula, w 0 (T), w 1 (T) are the weights of the foreground and the background respectively;

[0030] S213. Calculate the between-class variance σ 0 (T) that reflects the separation degree between different classes according to the foreground mean μ 1 (T) and the background mean μ 2 (T). The calculation formula is:

[0031] σ 2 (T) = w 0 (T) · w 1 (T) · (μ 0 (T) - μ 1 (T)) 2

[0032] S214. Traverse all possible thresholds T from 0 to 255 i , and calculate the between-class variance σ i under each threshold T 2 (T i );

[0033] S215. Select the threshold T 2 (T i ) that maximizes the between-class variance σ i as the global segmentation threshold T.

[0034] Further, in step S3, the specific process includes the following steps:

[0035] S31. Calculate the connected components of binary image bw1 and binary image bw2 and obtain the area S of each connected component i to form the dataset S;

[0036] S32. Take the median M of the dataset S and calculate the absolute difference D between the area S of each connected component i and the median M i to form the absolute deviation list D, where:

[0037] The expression of the median M is:

[0038] M = median(S)

[0039] The absolute difference D i is calculated by the formula:

[0040] D i = |S i - M

[0041] The expression of the absolute deviation list D is:

[0042] D = [D 1 , D 2 , …, D n , i ∈ n

[0043] In the formula, median(S) refers to the operation of taking the value at the middle position of the dataset S after sorting the dataset S from smallest to largest;

[0044] S33. Calculate the median of the absolute deviation list D to obtain the median absolute deviation MAD, that is:

[0045] MAD = median(D)

[0046] S34. Define the lower limit LB and upper limit UB for removing extreme values according to the median absolute deviation MAD, that is:

[0047] LB = M - MAD * k

[0048] UB = M + MAD * k

[0049] In the formula, k is a constant;

[0050] S35. Traverse the dataset S and remove the areas S that are less than the lower limit LB and greater than the upper limit UB i ;

[0051] S36. For the areas S of the remaining connected components after removing extreme values i′Calculate the new median, which is the particle reference area S O .

[0052] Further, in step S4, the specific process includes the following steps:

[0053] S41. Calculate the roundness C of the connected component i , and the calculation formula is:

[0054] C i = 4πS i / P 2

[0055] In the formula, S i represents the area of the i-th connected component; P represents the perimeter of the connected component;

[0056] S42. Filter the connected components according to the particle reference area S O and the roundness C i . When S O / S i > T S || C i < T C , then the connected component is a non-boundary abnormal area;

[0057] S43. Set all the pixel values included in the connected components corresponding to the non-boundary abnormal areas in the binary image bw1 and the binary image bw2 to 0;

[0058] S44. Perform an OR operation on the binary image bw1 and the binary image bw2 after removing the non-boundary abnormal areas to obtain a new binary image bw. The pixel bw(x, y) at the x-th column and the y-th row in the binary image bw is:

[0059]

[0060] In the formula, bw1(x, y) and bw2(x, y) are the pixels in the binary image bw1 and the binary image bw2 respectively.

[0061] Further, in step S5, the specific process includes the following steps:

[0062] S51. Generate a background image B with the same size as the particle image I based on the binary image bw, and set the initial gray value to the invalid value nan;

[0063] S52. Set the value of each pixel in the binary image bw to 0 to obtain a background image B with the same gray value as the particle image I, where the value of 0 represents the background and the value of 1 represents the particle;

[0064] The expression of the background image B is:

[0065]

[0066] Where bw(x,y) represents any pixel in the binary image bw;

[0067] S53, filling the remaining invalid values ​​nan in the background image B using B-spline interpolation according to the surrounding grayscale information.

[0068] Furthermore, in step S6, the specific process includes the following steps:

[0069] S61, calculate the difference image DI between the particle image I and the background image B, that is:

[0070] DI(x,y)=|I(x,y)-B(x,y)|

[0071] Where B(x,y) represents any pixel in the background image B;

[0072] S62, using a high-pass filter to remove low-frequency information in the difference image DI;

[0073] S63, using the difference threshold T d Identify the difference image DI that is greater than the difference threshold T d The pixel area is determined as the particle area, and the particle image I is obtained. * ;

[0074] Particle Image I * The expression is:

[0075]

[0076] S64, Particle Image I * Perform Gaussian filtering.

[0077] By means of the above technical solution, the present invention provides a particle image segmentation method for particle image velocimetry, which has at least the following beneficial effects:

[0078] 1. The method proposed in the present invention can remove interference while retaining the particle information in the image to the greatest extent. The interference factors in the particle image are identified and removed, and only valid particles are retained. This method has strong versatility and can provide more accurate and detailed data support for subsequent flow field research and analysis, and has broad application prospects.

[0079] 2. The present invention can solve and eliminate interference factors in the captured data, such as the influence of background objects, other light sources, etc. on PIV measurement, while removing interference and retaining the particle information in the image to the greatest extent, thereby identifying and removing interference factors in the particle image and retaining only valid particles. BRIEF DESCRIPTION OF THE DRAWINGS

[0080] The accompanying drawings described herein are used to provide a further understanding of the present application, and constitute a part of the present application. The schematic embodiments of the present application and their descriptions are used to explain the present application, and do not constitute an improper limitation of the present application. In the accompanying drawings:

[0081] Figure 1 is a schematic structural diagram of the two-dimensional PIV measurement device in the present invention;

[0082] Figure 2 is a schematic diagram of the particle image I collected in the present invention;

[0083] Figure 3 is a schematic diagram of the background image B generated from the particle image I in the present invention;

[0084] Figure 4 is a schematic diagram of the particle image I generated from the particle image I in the present invention * ;

[0085] Figure 5 is a schematic diagram of the particle image I with complex background collected in the present invention under a complex background scenario;

[0086] Figure 6 is a schematic diagram of the background image B generated from the particle image I with complex background in the present invention;

[0087] Figure 7 is a schematic diagram of the particle image I generated from the particle image I with complex background in the present invention * ; Detailed embodiments

[0088] In order to make the above objects, features, and advantages of the present invention more obvious and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. Thereby, a full understanding of how the present application uses technical means to solve technical problems and achieve technical effects can be obtained and implemented accordingly.

[0089] In the prior art, such as the Chinese invention patent with the publication number CN115170639A and the invention name of an image processing method for underwater three-dimensional particle image velocimetry of a robotic fish, it is proposed to preprocess the obtained tracer particle image, use the tracer particle image after two threshold segmentations to extract the bright part and the dark part of the fish body respectively, and then reconstruct the two-dimensional morphological characteristics of the entire robotic fish body, avoiding the incomplete extraction of the fish body caused by the different colors of the fish body and the single threshold setting.

[0090] The invention shows that the B, G, and R channels of the particle image are first segmented in turn, and the image with the highest contrast in each channel is selected. Then the particle image is median filtered, and the tracer particle image after median filtering is fused with the tracer particle image with the best contrast; the tracer particles are eliminated by reducing the gamma parameter. Based on the grayscale values ​​of the bright and dark parts of the fish body, a suitable threshold is given for thresholding and de-thresholding, and the binary images obtained after thresholding and de-thresholding are used as masks for the bright and dark parts of the fish body to achieve the extraction of the bright and dark parts of the fish body. Finally, the bright and dark parts are fused to determine the two-dimensional morphology of the entire fish body. However, the technical solution proposed by the invention has the following defects:

[0091] Firstly, for the particle images acquired by particle image velocimetry in the existing goldfish teaching experiment system, what is acquired is usually a grayscale image sequence, and it is impossible to obtain B, G, and R three-channel images.

[0092] Secondly, the binary image obtained after thresholding and de-thresholding is used as a mask for the bright and dark parts of the fish body. The premise is that the original image is globally thresholded. If the peak of the grayscale statistical distribution is not obvious, mis-segmentation is likely to occur, and the accurate edge of the fish body cannot be obtained.

[0093] In addition, this method is limited to a specific scenario (underwater robotic fish movement) and has poor scalability.

[0094] This embodiment aims to solve and eliminate the influence of interference factors (such as background objects, other light sources, etc.) in the captured data on PIV measurement, and to retain the particle information in the image to the greatest extent while removing the interference. Figures 1-7 This embodiment proposes a particle image segmentation method for particle image velocimetry, which can identify and remove interference factors in the particle image and retain only valid particles. This method is highly versatile and can provide more accurate and detailed data support for subsequent flow field research and analysis, and has broad application prospects. The method includes the following steps:

[0095] S1, obtain the particle image I located in the measurement area through the two-dimensional PIV measurement device; in this embodiment, the structural schematic diagram of the two-dimensional PIV measurement device is as follows Figure 1 As shown, it is composed of a high-speed camera, a surface laser, a water tank and other equipment. The water flow in the water tank is photographed with a high-speed camera to obtain high-quality image data. Specifically, when collecting particle images I, it is first necessary to arrange a suitable particle image velocity measurement environment through the above equipment, keep the water flow in the water tank flowing, and then adjust the surface laser settings to obtain the measurement plane, then spread the tracer particles to form the flow field environment to be measured, and finally collect the particle image I in the measurement area through the high-speed camera.

[0096] S2. Define the global segmentation threshold T and the local segmentation threshold T * , and segment the particle image I to obtain a binary image bw1 with the foreground and background separated, and a binary image bw2 of the particles in the particle image I that overlap with the background; in step S2, the specific process includes the following steps:

[0097] S21. Determine the optimal threshold T i as the global segmentation threshold T; in step S21, the specific process includes the following steps:

[0098] S211. Calculate the probability P(i) of each gray level i according to the gray level histogram of the particle image I, and the calculation formula is:

[0099]

[0100] In the formula, n i is the number of pixels of gray level i; N is the total number of pixels in the particle image I;

[0101] S212. Calculate the foreground mean μ 0 (T) and the background mean μ 1 (T) for the global segmentation threshold T. In this embodiment, the particle image I is divided into two categories, namely:

[0102] Foreground: All gray values are less than or equal to the global segmentation threshold T;

[0103] Background: All gray values are greater than the global segmentation threshold T.

[0104] Therefore, the calculation formulas for the foreground mean μ 0 (T) and the background mean μ 1 (T) are respectively:

[0105]

[0106] In the formula, w 0 (T), w 1 (T) are the weights of the foreground and the background respectively.

[0107] S213. Calculate the between-class variance σ 0 (T) that reflects the separation degree between different classes according to the foreground mean μ 1 (T) and the background mean μ 2 (T), and the calculation formula is:

[0108] σ 2 (T) = w 0 (T) · w 1 (T) · (μ 0 (T) - μ 1(T)) 2

[0109] S214. Traverse all possible thresholds T from 0 to 255 i , and calculate the between-class variance σ i under each threshold T 2 (T i );

[0110] S215. Select the threshold T that maximizes the between-class variance σ 2 (T i ) as the global segmentation threshold T. i

[0111] In this embodiment, by redefining the global segmentation threshold T and the local segmentation threshold T * , the target object in the image can be clearly distinguished from the background, and the contour and structure of the particle image can be highlighted, ensuring the accurate segmentation between the foreground and the background, as well as the overlapping particles in the background, thereby improving the accuracy of image processing. At the same time, it can process particle images under different lighting conditions and different contrasts, so as to obtain a more accurate and clearer binary image, providing strong support for subsequent steps such as contour extraction and generation of the background image B.

[0112] S22. Perform binary processing on the particle image I using the global segmentation threshold T to obtain a binary image bw1. In this embodiment, the particle image I is segmented using the global segmentation threshold T, and the purpose is to segment the foreground and the background in the particle image I. Therefore, the expression for the pixel bw1(x, y) in the x-th column and y-th row of the binary image bw1 is as follows:

[0113]

[0114] In the formula, I(x, y) represents any pixel point in the particle image I;

[0115] S23. Determine the optimal threshold T i * as the local segmentation threshold T * ;

[0116] S24. Perform binary processing on the particle image I using the local segmentation threshold T * to obtain a binary image bw2. The expression for the pixel bw2(x, y) in the x-th column and y-th row of the binary image bw2 is as follows:

[0117]

[0118] ​In this embodiment, the particle image I is binarized to segment the particles in the particle image I that overlap with the background. For the pixel I(x, y) in the x-th column and y-th row of the particle image I, the local segmentation threshold T is calculated according to the neighborhood pixels around this pixel (here, a sliding window with a default size of 51 is used) by the same principle as the method for determining the global segmentation threshold T above. * , and then the selected local segmentation threshold T * is used to convert the particle image I into a binary image bw2.

[0119] S3. Perform connected component calculation on the binary image bw1 and the binary image bw2, and remove extreme values to determine the particle reference area S O ; In this embodiment, the area S of each connected component is obtained by performing connected component calculation on the binary image bw1 and the binary image bw2 i . Since the particle sizes of the tracer particles are basically the same, the particle reference area S can be determined by removing extreme values O . The specific process includes the following steps:

[0120] S31. Perform connected component calculation on the binary image bw1 and the binary image bw2 and obtain the area S of each connected component i to form a data set S. Here, Python and the OpenCV library can be used to implement it. Use OpenCV to read the binary image bw1 and the binary image bw2, and then use the connectedComponentsWithStats function of OpenCV to find the connected components and obtain some statistical information of each connected component, including the area. Then extract the area S of each connected component from the returned statistical information i .

[0121] S32. Take the median M of the data set S and calculate the absolute difference D between the area S of each connected component i and the median M i to form an absolute deviation list D, where:

[0122] The expression of the median M is:

[0123] M = median(S)

[0124] The absolute difference D i is calculated by the formula:

[0125] D i = |S i - M

[0126] The expression of the absolute deviation list D is:

[0127] D = [D 1 , D2 , …, D n , i ∈ n

[0128] Where median(S) refers to the operation of taking the value at the middle position of the data set S after sorting S from smallest to largest;

[0129] S33. Calculate the median of the absolute deviation list D to obtain the median absolute deviation MAD, that is:

[0130] MAD = median(D)

[0131] S34. Define the lower limit LB and upper limit UB for removing extreme values according to the median absolute deviation MAD, that is:

[0132] LB = M - MAD * k

[0133] UB = M + MAD * k

[0134] Where k is a constant, and its value is 10 in this embodiment;

[0135] S35. Traverse the data set S, and remove the areas S that are less than the lower limit LB and greater than the upper limit UB i ;

[0136] S36. For the areas S of the remaining connected regions after removing extreme values i ′, calculate a new median, and this median is the particle reference area S O .

[0137] S4. Use the particle reference area S O to filter the connected regions that are non-boundary abnormal regions to obtain a new binary image bw; in step S4, the specific process includes the following steps:

[0138] S41. Calculate the roundness C of the connected region i , and the calculation formula is:

[0139] C i = 4πS i / P 2

[0140] Where S i represents the area of the i-th connected region; P represents the perimeter of the connected region;

[0141] S42. Filter the connected regions according to the particle reference area S O and the roundness C i . When S O / S i > T S || C i < T CIf so, the connected region is a non-boundary abnormal region. In this embodiment, the value of T S is 10, and the value of T C is 0.5.

[0142] S43. Set all the pixel values included in the connected regions corresponding to the non-boundary abnormal regions in the binary image bw1 and the binary image bw2 to 0;

[0143] S44. Perform an OR operation on the binary image bw1 and the binary image bw2 after removing the non-boundary abnormal regions to obtain a new binary image bw. At this time, the effective region values in the binary image bw1 and the binary image bw2 are 1, the invalid region and the non-boundary abnormal region values are 0. For the pixel bw(x, y) in the x-th column and y-th row of the binary image bw:

[0144]

[0145] In the formula, bw1(x, y) and bw2(x, y) are the pixels in the binary image bw1 and the binary image bw2 respectively.

[0146] S5. Generate a background image B that is consistent with the gray value of the particle image I based on the binary image bw; In step S5, the specific process includes the following steps:

[0147] S51. Generate a background image B with the same size as the particle image I based on the binary image bw, and set the initial gray value to the invalid value nan;

[0148] S52. Set the value of each pixel in the binary image bw to 0 to obtain a background image B that is consistent with the gray value of the particle image I, where the value of 0 represents the background and the value of 1 represents the particle;

[0149] The expression of the background image B is:

[0150]

[0151] In the formula, bw(x, y) represents any pixel in the binary image bw.

[0152] S53. Fill the remaining invalid values nan in the background image B according to the surrounding gray information using B-spline interpolation. More specifically, B-spline interpolation is often used for data interpolation and filling missing values. Its interpolation model is a piecewise polynomial function, which consists of several local polynomials that are connected at given nodes. When constructing the B-spline interpolation model, it is necessary to set the nodes (key points defined by the B-spline curve) and the degree (degree of the local polynomial). The nodes and the degree affect the smoothness and accuracy of the interpolation. In this embodiment, the default uniform nodes are used, and the degree is 3.

[0153] S6. Segment the particle image I according to the background image B to obtain the particle image I * ; In step S6, the specific process includes the following steps:

[0154] S61. Calculate the difference image DI between the particle image I and the background image B, that is:

[0155] DI(x,y) = |I(x,y) - B(x,y)|

[0156] In the formula, B(x,y) represents any pixel in the background image B;

[0157] S62. Use a high-pass filter to remove the low-frequency information in the difference image DI and improve the clarity of the particle signal;

[0158] S63. Set an appropriate difference threshold T d (generally 10 - 20, the default value is 20), identify the pixel regions in the difference image DI that are greater than the difference threshold T d and determine them as particle regions, and obtain the particle image I accordingly * ;

[0159] Particle image I * The expression is:

[0160]

[0161] S64. Perform Gaussian filtering on the particle image I * to ensure a smooth transition between the particle edges and the background.

[0162] The method proposed in this embodiment can, in PIV measurement, remove interference while maximizing the retention of particle information in the image. Identify and remove the interference factors in the particle image, and only retain the effective particles. This method has strong versatility and can provide more accurate and detailed data support for the subsequent research and analysis of the flow field, and has broad application prospects.

[0163] Those of ordinary skill in the art can understand that all or part of the steps in implementing the method of the above embodiment can be completed by instructing relevant hardware through a program. Therefore, this application can adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, this application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk memories, CD-ROMs, optical memories, etc.) containing computer-usable program codes.

[0164] The above embodiments have introduced the present invention in detail. Specific examples are used in this article to elaborate on the principles and embodiments of the present invention. The description of the above embodiments is only used to help understand the method and its core idea of the present invention; at the same time, for those of ordinary skill in the art, according to the idea of the present invention, there will be changes in the specific implementation manners and application scopes. In summary, the content of this specification should not be construed as a limitation to the present invention.

Claims

1. A particle image segmentation method for particle image velocimetry, characterized in that: The method comprises the following steps: S1, obtaining a particle image I located in a measurement area through a two-dimensional PIV measurement device; S2. Define the global segmentation threshold T and the local segmentation threshold T * , and segment the particle image I to obtain a binary image bw1 in which the foreground and the background are separated, and a binary image bw2 of the particles in the particle image I that overlap with the background is segmented; S3, calculate the connected domain of the binary image bw1 and the binary image bw2, and remove the extreme values ​​to determine the particle reference area S O ; S4, using particle reference area S O Filter the connected domains that are not boundary abnormal areas to obtain a new binary image bw; S5, generating a background image B consistent with the grayscale value of the particle image I based on the binary image bw; S6. Segment the particle image I according to the background image B to obtain the particle image I * .

2. The particle image segmentation method according to claim 1, characterized in that: The process of obtaining the particle image I comprises: A particle image velocimetry environment was built using a high-speed camera, a surface laser, and a water tank. The water flow in the water tank was kept flowing, and then the surface laser settings were adjusted to obtain the measurement plane. Subsequently, tracer particles were spread to form a flow field environment to be measured. Finally, a high-speed camera was used to collect the particle image I in the measurement area.

3. The particle image segmentation method according to claim 1, characterized in that: In step S2, the specific process includes the following steps: S21. Determine the optimal threshold T i As the global segmentation threshold T; S22, using the global segmentation threshold T to perform binarization processing on the particle image I to obtain a binary image bw1, the expression of the pixel bw1(x, y) in the xth column and yth row in the binary image bw1 is: S23. Determine the optimal threshold T i * As the local segmentation threshold T * ; S24, using local segmentation threshold T * The particle image I is binarized to obtain a binary image bw2. The expression of the pixel bw2(x,y) in the xth column and yth row in the binary image bw2 is: Where I(x,y) represents any pixel in the particle image I.

4. The particle image segmentation method according to claim 3, characterized in that: In step S21, the specific process includes the following steps: S211, calculate the probability P(i) of each gray level i according to the gray histogram of the particle image I, and the calculation formula is: Where n i is the number of pixels at gray level i; N is the total number of pixels in particle image I; S212. Calculate the foreground mean μ0(T) and background mean μ1(T) at different gray levels i according to the probability P(i). The calculation formulas for the foreground mean μ0(T) and the background mean μ1(T) are: In the formula, w0(T) and w1(T) are the weights of the foreground and background respectively; S213, calculate the inter-class variance σ to reflect the degree of separation between different classes based on the foreground mean μ0(T) and the background mean μ1(T) 2 (T), the calculation formula is: σ 2 (T)=w0(T)·w1(T)·(μ0(T)-μ1(T)) 2 S214, traverse all possible threshold values ​​T from 0 to 255 i , and calculate each threshold T i The between-class variance σ 2 (T i ); S215, select the between-class variance σ 2 (T i )The maximum threshold T i As the global segmentation threshold T.

5. The particle image segmentation method according to claim 1, characterized in that: In step S3, the specific process includes the following steps: S31, calculate the connected domains of the binary image bw1 and the binary image bw2 and obtain the area S of each connected domain i Construct a data set S; S32. Take the median M of the data set S, and calculate the area S of each connected domain based on the median M. i The absolute difference D from the median M i Construct a list of absolute deviations D, where: The expression of median M is: M=median(S) Absolute difference D i The calculation formula is: D i =|S i -M| The expression of the absolute deviation list D is: D=[D1,D2,…,D n ],i∈n In the formula, median(S) refers to the operation of taking the middle value of the data set S after sorting the data set S from small to large; S33, calculate the median of the absolute deviation list D, and obtain the median absolute deviation MAD, that is: MAD=median(D) S34. Define the lower limit LB and upper limit UB for removing extreme values ​​according to the median absolute deviation MAD, that is: LB=M-MAD*k UB=M+MAD*k Where k is a constant; S35, traverse the data set S, and find the area S that is smaller than the lower limit LB and larger than the upper limit UB i Remove; S36, the area S of the connected domain remaining after removing the extreme values i ' Calculate the new median, which is the particle reference area S O .

6. The particle image segmentation method according to claim 1, characterized in that: In step S4, the specific process includes the following steps: S41. Calculate the circularity C of the connected domain i , the calculation formula is: C i =4πS i / P 2 In the formula, S i represents the area of ​​the i-th connected domain; P represents the perimeter of the connected domain; S42, according to the particle reference area S O and roundness C i Filter the connected domain. When S O / S i >T S ||C i <T C When , the connected domain is a non-boundary abnormal area; S43, setting all pixel values ​​contained in the connected domains corresponding to the non-boundary abnormal regions in the binary image bw1 and the binary image bw2 to 0; S44, the binary image bw1 and the binary image bw2 after removing the non-boundary abnormal area are combined to obtain a new binary image bw, and the pixel bw(x, y) in the xth column and yth row of the binary image bw is: Where bw1(x,y) and bw2(x,y) are pixels in the binary image bw1 and the binary image bw2 respectively.

7. The particle image segmentation method according to claim 1, characterized in that: In step S5, the specific process includes the following steps: S51, generating a background image B of the same size as the particle image I based on the binary image bw, and setting the initial grayscale value to an invalid value nan; S52, setting the value of each pixel in the binary image bw to 0, and obtaining a background image B having the same grayscale value as the particle image I, wherein a value of 0 represents the background, and a value of 1 represents the particle; The expression of the background image B is: Where bw(x,y) represents any pixel in the binary image bw; S53, filling the remaining invalid values ​​nan in the background image B using B-spline interpolation according to the surrounding grayscale information.

8. The particle image segmentation method according to claim 1, characterized in that: In step S6, the specific process includes the following steps: S61, calculate the difference image DI between the particle image I and the background image B, that is: DI(x,y)=|I(x,y)-B(x,y)| Where B(x,y) represents any pixel in the background image B; S62, using a high-pass filter to remove low-frequency information in the difference image DI; S63, using the difference threshold T d Identify the difference image DI that is greater than the difference threshold T d The pixel area is determined as the particle area, and the particle image I is obtained. * ; Particle Image I * The expression is: S64, Particle Image I * Perform Gaussian filtering.

Images