Erosion pit identification and parameter extraction method based on three-dimensional extreme point

Abstract

The invention discloses a corrosion pit identification and parameter extraction method based on a three-dimensional extreme point, and belongs to the field of corrosion detection of metal structures, and the method comprises the following steps: S1, obtaining surface point cloud plane original data of a corrosion test piece through three-dimensional laser scanning, obtaining standard two-dimensional height grid curved surface data through point cloud plane coordinate system conversion and data repair based on the original data; s2, extracting effective extreme points of the curved surface by selecting a neighborhood detection range, extreme point saliency and extreme point spacing; s3, Gaussian smoothing processing is performed on the grid curved surface data, and a curved surface global gradient is calculated through a Sobel operator based on the smoothed curved surface; and S4, obtaining an extreme point boundary based on the segmented label matrix, and carrying out corresponding statistical processing on region data to which each extreme point belongs. By means of the method, automatic identification and three-dimensional parameter accurate measurement of the corrosion area can be achieved, and a new solution is provided for metal structure corrosion evaluation.

Description

TECHNICAL FIELD

[0001] The present application relates to a kind of based on three-dimensional extreme point's etch pit identification and parameter extraction method. BACKGROUND

[0002] In industrial infrastructure operation and maintenance, major project safety monitoring and other fields, corrosion detection of metal structure is the key link to ensure the safe operation of equipment. Corrosion not only causes material mechanical properties degradation, structure load capacity decline, but also can cause sudden fracture, leakage and other safety accidents, causing huge economic losses and environmental risks. Current corrosion detection mainly relies on manual visual inspection, two-dimensional image analysis and other traditional methods, which have low detection accuracy, strong subjectivity and cannot quantify the three-dimensional topography parameters (such as depth, volume, surface area) of corrosion pits. Especially for local corrosion pits on the surface of complex components, traditional technology cannot achieve high-precision positioning and multi-dimensional feature extraction, resulting in corrosion degree evaluation deviation, which directly affects the accuracy of structure residual life prediction. Therefore, developing detection technology that can accurately identify and quantify the three-dimensional features of corroded surface has important practical significance for ensuring the safe operation of major engineering structures and reducing maintenance costs. SUMMARY

[0003] The technical problem to be solved by the present application is to provide a three-dimensional extreme point-based etch pit identification and parameter extraction method to realize automatic identification of corrosion area and accurate measurement of three-dimensional parameters, and provide a new solution for metal structure corrosion evaluation.

[0004] The technical solution adopted by the present application to solve the technical problem is: a three-dimensional extreme point-based etch pit identification and parameter extraction method, comprising the following steps: S1: obtaining the original data of corrosion specimen surface point cloud by three-dimensional laser scanning, and obtaining standard height grid surface data by point cloud plane coordinate system conversion and data repair based on the original data; S2: extracting effective extreme points from the surface based on the standard height grid surface data by selecting neighborhood detection range, extreme point saliency and extreme point spacing; S3: performing Gaussian smoothing on the grid surface data and calculating the global gradient of the surface by Sobel operator based on the smoothed surface, determining the background area by distance transformation, and obtaining the region boundary by watershed algorithm based on the corrected gradient with the extreme points obtained in S2 as local minimum gradient; S4: obtaining the extreme point boundary based on the segmented label matrix, performing corresponding statistical processing on the region data of each extreme point, extracting the etch pit depth and diameter corresponding to the extreme point region, and counting the overall etch pit data distribution.

[0005] Further, in step S1, the surface point cloud original data of the corrosion specimen is obtained by three-dimensional laser scanning, the original data points are fitted to a plane equation to solve a tilt plane, a rotation matrix is constructed according to the normal vector of the fitting plane to rotate the tilt plane for the first time, the main direction of the plane is determined by principal component analysis to construct a rotation matrix to rotate the plane for the second time to align the x and y coordinate axes, the regular grid data parallel to the x and y coordinate axes is generated by bilinear interpolation from the discrete plane coordinate data after rotation, and the original data missing values are filled in the nearest neighbor manner, so that the regular standard grid data is finally obtained, and the specific operation steps are as follows: (1) Original point cloud acquisition: the surface point cloud original data of the corrosion specimen is obtained by non-contact three-dimensional laser scanning , and an inclined original plane is obtained. (2) Inclined plane equation fitting: assuming that the plane equation constructed by the original data points is , the original coordinate data of the corrosion surface is obtained by non-contact three-dimensional laser scanning , and the original surface matrix is . According to the least square method, the plane coefficients can be solved: ; wherein is the original height vector. (3) Inclined plane rotation: according to the normal vector of the original plane equation and the normal vector of the x and y plane, the rotation axis and the rotation angle of the two normal vectors are calculated, and the rotation matrix is obtained by the rotation axis and the rotation angle , so as to correct the tilt of the plane. , ; , ; , wherein , ; the point coordinates after rotation , wherein is the original data point obtained by scanning, and the plane is translated in the Z direction so that the lowest point of the plane is , and the initial tilt original data is converted to the horizontal x-y plane coordinate system. (4) Second plane rotation: the horizontal plane point cloud data after the first rotation and translation is extracted as two-dimensional coordinates ; compute the covariance matrix of two-dimensional coordinates , eigenvalues of the covariance matrix are obtained by eigenvalue decomposition , and the corresponding eigenvectors , , a rotation matrix can be constructed ; ; , ; ; where , is the mean of the x, y coordinates of the point cloud in the plane, n is the number of points, is the main direction representing the longest direction of data extension, represents the direction perpendicular to the main direction; After rotating the two-dimensional coordinates in the x-y plane , the coordinates are obtained, and the center of the point cloud in the plane is aligned with the origin, the main direction of the point cloud is aligned with the x-axis, and the secondary direction of the point cloud is aligned with the y-axis to realize the regular distribution of the data in the plane; (5) Grid data fitting and missing value filling: according to the x, y axis range of the point cloud coordinates after the second rotation and the grid spacing, the number of grid nodes in the x, y direction is determined: , ; , ; , ; , ; where , is the grid accuracy input actively; , is the number of grid nodes in the x, y direction; Since the point cloud data is discrete, the height value of each grid node will be calculated by bilinear interpolation, based on the weighted average of the four adjacent point clouds around the node, where the upper left corner point (x1, y1, z1), the upper right corner point (x3, y1, z3), the lower left corner point (x1, y3, z3), and the lower right corner point (x3, y3, z4): x-direction interpolation (fixed y): ; ; y direction interpolation (fixed x): ; Since the initial scanning data may be non-rectangular or the boundaries are uneven, the initial grid boundary obtained may not be complete, and there will be missing values in the rectangular grid matrix. The missing values are recorded as and filled in the nearest neighbor way. Specifically, first, fill the missing values along the column direction of the grid data: for each column , traverse the row index , if , search for the nearest non-missing value in the column , fill with ; then fill it again along the row direction in the above manner; finally, obtain the standard curved surface grid vertex height data .

[0006] Further, in step S2, the local neighborhood range is defined by analyzing the control parameters, the candidate extreme value points are screened based on the neighborhood extreme value, mean difference and saliency, the spatial distance constraint is removed, the final extreme value point coordinate set of the grid surface data is generated and used for subsequent feature analysis, and the specific steps are as follows: (1) Define control parameters: minimum saliency (minimum difference value between extreme value point and neighborhood), minimum distance (minimum spatial distance between extreme value points), neighborhood size (defining the size of local range row x column); (2) Define local neighborhood: let the three-dimensional curved surface data surface_data be a matrix , calculate the neighborhood radius for each detection point : , , define the neighborhood range: ; ; Where R represents a real number, M and N are the row and column numbers corresponding to the grid data; (3) Screen extreme value points: record the neighborhood sub-matrix , calculate the neighborhood mean value and saliency for each non-boundary point: ; ; If it meets: ; , mark For candidate maximum value points, corresponding coordinates are stored in a set , corresponding height values are stored in ; (4) Maximum value point constraint de-duplication: spatial de-duplication is performed on the maximum value point candidates, and is sorted in descending order, and the sorted points are calculated with the Euclidean distance of the reserved points : ; Where the reserved points are initially empty, and when the candidate points are calculated one by one, if , then is added to ; otherwise, it is removed, and the final maximum value point set is output after de-duplication .

[0007] Further, in step S3, after smoothing the grid surface data by Gaussian filtering, the gradient amplitude of the grid data is calculated based on the Sobel operator, the background area and the foreground area of the gradient graph are determined by expanding the data gradient amplitude through distance transformation, and the extreme value points are forcibly marked as local minima of the gradient graph. The label matrix is generated by applying the watershed algorithm based on the modified gradient graph to realize region segmentation, and the image segmentation of the background and foreground areas is completed. The specific steps are as follows: (1) Data smoothing processing: in order to reduce the interference of noise on gradient calculation, two-dimensional Gaussian filtering is used to smooth the original surface data, and the smoothed data is the convolution of the original surface data and the Gaussian kernel: ; The Gaussian kernel function is: , where , is the offset of the pixel in the kernel relative to the center; (2) Surface gradient calculation: the gradient amplitude of the smoothed surface data is calculated by using the Sobel operator, and each gradient component can be obtained by convolution calculation: , ; Where the Sobel operator kernel is defined as , ; According to the gradient components, the gradient amplitude can be calculated as: ; (3) Background calculation: first, mark the spatial position of the extreme value points, define the mark matrix ​(Same size as the original mesh surface data), initialized as an all-zero matrix; for each extreme point... ,set up The "background region" and the "foreground region" containing the extreme points are determined through distance transformation; the Euclidean distance from each pixel to the nearest extreme point is calculated. ; in, Let the coordinates be the k-th extreme point. When the pixel distance is greater than the threshold multiplied by the maximum distance, that is... Mark as background Update the marker simultaneously , where yuzhi is the input parameter, which controls the size of the background area and is usually set to 0.4; (4) Gradient magnitude correction: In order to constrain the "watering start point" and "region expansion range" of the watershed algorithm, and to ensure that the watershed algorithm "floods" the region from the extreme point, while avoiding the over-segmentation problem of the unconstrained watershed algorithm, the foreground extreme point and the background region will be forced to be set as local minima of the gradient topographic map through mathematical constraints, so that both become watering points at the same time. Treating the gradient magnitude as terrain height, we have: ; From the tag matrix Extract the coordinate set of all foreground extreme points: ; Where K is the total number of foreground extreme points; For each extreme point Define its domain as and all pixels in the neighborhood Perform gradient correction: ; in It is a very small positive number to ensure that the gradient value of the extreme point itself is strictly less than that of other points in the neighborhood; This represents all pixels in the neighborhood except for the extreme points; Similarly, the background region is corrected by the gradient map and forced into a local minimum region of gradient terrain, serving as the "background water injection starting point" of the watershed algorithm; From the label matrix, we can see that This is the background region, denoted as B. A region-wide gradient reduction strategy is applied to background region B to correct the gradient: ; in This is the original gradient map. Represents the minimum gradient value in the entire image. It is a very small positive number; The corrected background region becomes the flattest region globally, which can be regarded as the sea level in the watershed algorithm; The above steps complete the correction of the gradient assignment of the original data. (5) Watershed boundary division: Corrected gradient map The initial set of local minima is the union of the foreground and background regions: ; in (The lowest point within the background region ensures that the background is a single connected minimum region.) The initial region and label matrix are now initialized for the subsequent watershed algorithm: Foreground area: for each initial region (Includes only the extreme points themselves); Background area: Initial area (The entire background area is used as the initial area); Tag matrix: ; in Ensure that the background label and the foreground label are not duplicated; Flood fill the corrected gradient graph based on a priority queue: 1. Water level parameters and event queue: Water level range: ,in (The initial water level is the gradient value of the background area) (Maximum gradient); Event queue: A queue Q is used to store pixels to be processed, arranged in ascending order of gradient value h(x,y). The simulation simulates water level rising from low to high. Initially, all minimum points are added to the queue. ; 2. Expansion of the flood area: Pixels are processed in ascending order of gradient values ​​(simulating water levels rising from low to high). For each pixel (x, y), a priority queue is used. Extract and process: First, extract the pixel (x,y) with the smallest h(x,y) from Q, and let the current water level t = h(x,y); then, start neighborhood expansion based on the minimum water level, and expand the eight neighborhoods of (x,y). Pixels Perform marking and determination: If this is not marked Then It is assigned to the current region and added to the queue to await processing of its neighbors: ; If marked and label is different , mark as pixel boundary: ; When the priority queue Q is empty, all pixels are processed, at this time: foreground region: ; background region: ; label matrix: ; wherein represents the other regions of the foreground region except the boundary points, represents the other regions of the background region except the boundary points.

[0008] Further, in the step S4, based on the label matrix obtained by the watershed algorithm, first, the boundary coordinate set corresponding to each extreme point is obtained; then the depth and diameter of the etch pit are extracted from each extreme point region through statistical processing, and the etch pit depth characteristics of the overall etching morphology are obtained, and the specific steps are as follows: (1) extreme point label matrix boundary extraction: according to the label matrix obtained in S3, the boundary pixel coordinates of the label region corresponding to each extreme point are extracted, and the three-dimensional boundary coordinates are generated by combining the original grid data height value. For the first region (label is ), the contour pixel coordinates are extracted: ; (2) etch pit feature extraction: according to the contour pixel coordinates combined with the original curved surface height data , the etch pit boundary coordinates can be generated: ; According to the extreme point and the boundary coordinates , the etch pit depth and the etch pit diameter can be determined. The maximum and minimum values of the etch pit boundary x and y directions are recorded as x kmin , x kmax , y kmin , y kmax , then ; , and the corresponding etch pit depth-diameter ratio can be expressed as ; (3) etch pit distribution statistics: interval distribution probability statistics is performed on the etch pit depth and depth-diameter ratio obtained in (2), and the probability distribution density histogram is obtained; then the data is fitted to a normal distribution, and the corresponding data distribution is obtained.​

[0009] The beneficial effects of the present application are that the irregular point cloud data obtained by three-dimensional scanning can be converted into regular grid height matrix data, and the data can be effectively screened and the corrosion pit detection under different corrosion stages can be effectively adapted by the extreme point multi-parameter control and the watershed boundary segmentation method based on the extreme point, thereby providing a quantitative basis for corrosion evaluation and helping researchers to conduct corrosion evolution research and life prediction. BRIEF DESCRIPTION OF DRAWINGS

[0010] Figure 1 is a three-dimensional extreme point-based pit recognition and parameter extraction method guide map; Figure 2 is a plane view of the original data point cloud; Figure 3 is a first rotation plane view of the inclined plane; Figure 4 is a second rotation main direction alignment view of the point cloud plane; Figure 5 is a point cloud plane grid view with data missing; Figure 6 is a point cloud plane grid view after missing value completion; Figure 7 is a three-dimensional extreme point extraction schematic view; Figure 8 is a raw data grayscale view and a filtered data grayscale view; Figure 9 is a data foreground and background division relying on the extreme point; Figure 10 is a corrected gradient grayscale view; Figure 11 is a label matrix obtained after division by the watershed algorithm; Figure 12 is a three-dimensional extreme point pit boundary division schematic view obtained according to the label matrix; Figure 13 is a pit depth-diameter ratio and pit depth data histogram and distribution. DETAILED DESCRIPTION

[0011] The present application will be further described below in combination with the drawings and examples.

[0012] With reference to Figure 1 , the present embodiment provides a three-dimensional extreme point-based pit recognition and parameter extraction method, including the following steps: S1: original data acquisition and preprocessing (1) Original point cloud acquisition: taking a Q355 steel plate corroded for 30 days under neutral salt spray as an example, the original data of the point cloud on the surface of the corrosion test piece is obtained by non-contact three-dimensional laser scanning , the obtained inclined original plane is as shown in Figure 2 ; (2) Inclined plane equation fitting: assuming that the plane equation constructed by the original data points is , the original coordinate data of the corroded surface obtained by the non-contact three-dimensional laser scanning is , the original surface matrix is , and the plane coefficients can be solved according to the least square method: ; wherein, is the original height vector; (3) Inclined plane rotation: according to the normal vector of the original plane equation and the normal vector of the x, y plane, the rotation axis and the rotation angle of the two normal vectors are calculated, and the rotation matrix can be obtained through the rotation axis and the rotation angle , wherein , ; the point coordinates after rotation are , wherein is the original data point obtained by scanning, and the plane is translated in the Z direction so that the lowest point of the plane is , and the initial inclined plane correction is completed, as shown in Figure 3 ; (4) Second plane rotation: the horizontal plane point cloud data after the first rotation and translation is extracted as two-dimensional coordinates ; the covariance matrix of the two-dimensional coordinates is calculated , wherein , are the mean values of the x, y coordinates of the point cloud in the plane, and n is the number of points; the eigenvalues and the corresponding eigenvectors , are obtained by eigenvalue decomposition of the covariance matrix, wherein is the main direction representing the longest direction of the data extension, represents the direction perpendicular to the main direction, the main direction is taken as the x-axis target direction, and the rotation matrix is constructed; the two-dimensional coordinates in the x-y plane are rotated to obtain the coordinates , so that the center of the plane point cloud is aligned with the origin, the main direction of the point cloud is aligned with the x axis, and the secondary direction of the point cloud is aligned with the y axis, so that the regularized distribution of the data in the plane is realized, as shown in Figure 4 . (5) Grid data fitting and missing value imputation: The number of grids in the x and y directions is determined based on the x and y axis ranges of the point cloud coordinates after the second rotation and the grid spacing. , ; , ; , ; , ; in , The mesh accuracy for active input is 0.5mm; , The number of grid nodes in the x and y directions; Because point cloud data is discrete, each grid node height value The height will be calculated using bilinear interpolation, based on the weighted average of the four neighboring point clouds surrounding the node, where the top left point is (x1, y1, z1), the top right point is (x3, y1, z3), the bottom left point is (x1, y3, z3), and the bottom right point is (x3, y3, z4). x-direction interpolation (with fixed y): ; ; y-direction interpolation (fixed x): ; Because the initial scan data may not be rectangular or may have irregular boundaries, the obtained initial mesh boundary may not be complete. Missing values ​​will exist in the rectangular mesh matrix; these missing values ​​will be denoted as... And fill according to the nearest neighbor method, specifically: first, fill missing values ​​in the grid data along the column direction: for each column traverse row indexes ,like Then search for the distance in that column. The most recent non-missing value ,Will Fill as Then, perform a second filling along the row direction as described above; finally, obtain the vertex height data of the standard curved surface mesh. like Figure 6 As shown; S2: Extreme Point Extraction and Filtering (1) Define control parameters: Set the specific set of control parameters Minimum significance (Controlling the difference between the extreme value point and the neighborhood), minimum distance (Controlling the minimum space distance between the extreme value points), neighborhood size (Defining the size of the local range row x column); (2) Define the local neighborhood: Let the three-dimensional surface data surface_data be a matrix , calculate the neighborhood radius for each detection point , , define the neighborhood range: ; ; where R represents a real number, M and N are the row and column numbers corresponding to the grid data; (3) Screening maximum value points: Let the neighborhood sub-matrix , calculate the neighborhood mean value for each non-boundary point and the saliency , if: ; , mark as a candidate maximum value point, and store the corresponding coordinates in the set , and store the corresponding height value in ; (4) Maximum value point constraint de-duplication: Spatial de-duplication is performed on the maximum value point candidate value, and is sorted in descending order, and for the sorted point , the Euclidean distance between the point and the retained point is calculated: ; where the retained point is initially an empty set, and when the candidate point is calculated Euclidean distance one by one, if , then is added to ; otherwise, it is removed, and the final maximum value point set is output after de-duplication , and the specific maximum value point position display is shown in Figure 7 ; S3: Watershed image segmentation (1) Data smoothing processing: In order to reduce the interference of noise on gradient calculation and to preserve the true trend of the data, a two-dimensional Gaussian filter is used to smooth the original surface data, and the Gaussian kernel function is: where , is the offset of the pixel in the kernel relative to the center; the smoothed data is the convolution of the original surface data and the Gaussian kernel: , the gray scale diagrams before and after correction are as follows​​Figure 8 As shown; (2) Surface gradient calculation: The smoothed surface data is calculated using the Sobel operator. The gradient magnitude; where the Sobel operator kernel is defined as , The gradient components can be obtained through convolution: , ; The gradient magnitude can be determined based on the gradient components. calculate: ; (3) Background calculation: First, mark the spatial location of the extreme points and define the marking moments. (Same size as the original mesh surface data), initialized as an all-zero matrix; for each extreme point... ,set up The "background region" and the "foreground region" containing the extreme points are determined through distance transformation; the Euclidean distance from each pixel to the nearest extreme point is calculated. ; in( (The coordinates of the k-th extreme point) When the pixel distance is greater than the threshold × the maximum distance, that is... Mark as background Update the marker simultaneously Where yuzhi is the input parameter, controlling the size of the background area, usually set to 0.4, as shown in the following example. Figure 9 As shown, The regions in the matrix where the elements are not zero are displayed in white, including extreme points and the background; The regions in the matrix with elements of 0 are displayed in black, representing the neighborhood of the extreme point; (4) Gradient magnitude correction: In order to constrain the "watering start point" and "region expansion range" of the watershed algorithm, and to ensure that the watershed algorithm "floods" the region from the extreme point, while avoiding the over-segmentation problem of the unconstrained watershed algorithm, the foreground extreme point and the background region will be forced to be set as local minima of the gradient topographic map through mathematical constraints, so that both become watering points at the same time. Treating the gradient magnitude as terrain height, we have: ; From the tag matrix Extract the coordinate set of all foreground extreme points: ; Where K is the total number of foreground extreme points; For each extreme point Define its domain as All pixels in the neighborhood Perform gradient correction: ; Where is a very small positive number to ensure that the gradient value of the extreme point itself is strictly less than other points in the neighborhood; represents other pixels in the neighborhood except the extreme point; Similarly, the background region is forced to be a local minimum region of the gradient terrain through gradient correction, as the "background water injection starting point" of the watershed algorithm; Through the marking matrix That is, the background region, now the background region is denoted as B, and the gradient of the background region B is corrected using the overall region depression strategy: ; Where is the original gradient map, represents the minimum gradient value in the entire image, is a very small positive number; The corrected background region becomes the global lowest flat region, which can be regarded as the sea level in the watershed algorithm; Through the above work, the correction of the original data gradient assignment is completed, and the corrected gradient gray image is shown in Figure 10 ; (5) Watershed boundary segmentation: The initial minimum point set of the corrected gradient map is the union of foreground and background regions: ; Where (the lowest value point in the background region, to ensure that the background is a single connected minimum value region); Now initialize the initial region and label matrix to proceed with the following watershed algorithm: Foreground region: for each initial region (only contains the extreme point itself); Background region: initial region (the entire background region as the initial region); Label matrix: ; Where , to ensure that the background label and the foreground label are not repeated; Based on the priority queue, flood fill the corrected gradient map: 1. Water level parameter and event queue: Water level interval: , where (The initial water level is the gradient value of the background area) (Maximum gradient); Event queue: A queue Q is used to store pixels to be processed, arranged in ascending order of gradient value h(x,y). The simulation simulates water level rising from low to high. Initially, all minimum points are added to the queue. ; 2. Expansion of the flood area: Pixels are processed in ascending order of gradient values ​​(simulating water levels rising from low to high). For each pixel (x, y), a priority queue is used. Extract and process: First, extract the pixel (x,y) with the smallest h(x,y) from Q, and let the current water level t = h(x,y); then, start neighborhood expansion based on the minimum water level, and expand the eight neighborhoods of (x,y). Pixels Perform marking and determination: If this is not marked Then It is assigned to the current region and added to the queue to await processing of its neighbors: ; If already marked and the labels are different Then it is marked as a pixel boundary: ; When priority queue Q is empty, all pixels are processed. At this time: Foreground area: ; Background area: ; Tag matrix: ; in This represents the foreground region excluding the boundary points. This represents the background area excluding the boundary points. S4: Calculation of pitting parameters (1) Boundary extraction of the extreme point label matrix: Based on the label matrix obtained in S3, the boundary pixel coordinates of the label region corresponding to each extreme point are extracted and combined with the height value of the original grid data to generate three-dimensional boundary coordinates. Areas (labeled as) Extract its contour pixel coordinates. The boundaries of each corrosion pit are as follows Figure 12 The red lines shown represent the boundaries of each extreme point: ; (2) Pits feature extraction: Based on the contour pixel coordinates combined with the original surface height data This will generate the coordinates of the erosion pit boundary: ; Based on extreme points and boundary coordinates The depth of the erosion pit can then be determined. Diameter of the erosion pit Let the maximum and minimum values ​​of the pit boundary in the x and y directions be denoted as x, y, y, respectively. kmin x kmax y kmin y kmax So there are ; The corresponding pit depth-to-diameter ratio can be expressed as: ; (3) Distribution statistics of pits: The pit depth and depth-to-diameter ratio obtained in (2) are statistically analyzed for interval distribution to obtain a probability distribution density histogram; then, the data are fitted with a normal distribution to obtain the following results. Figure 13 The corresponding data distribution is shown below. Figure 13 The probability distribution density shown on the vertical axis is the interval probability of the data divided by the unit interval length.

[0013] Although the preferred embodiments of the present invention have been described above in conjunction with the accompanying drawings, the present invention is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make many other modifications under the guidance of the present invention without departing from the spirit and scope of the claims, and all of these modifications are within the scope of protection of the present invention.

Claims

1. A method for identifying erosion pits and extracting parameters based on three-dimensional extreme points, characterized in that, Includes the following steps: S1: Obtain the original point cloud data of the corrosion specimen surface through three-dimensional laser scanning, and obtain standard height grid surface data based on the original data by point cloud plane coordinate system transformation and data repair. S2: Based on standard height grid surface data, effective extreme points are extracted from the surface by selecting the neighborhood detection range, extreme point significance, and extreme point spacing; S3: Gaussian smoothing is applied to the mesh surface data, and the global gradient of the surface is calculated using the Sobel operator based on the smoothed surface. The background region is determined by distance transformation. Based on the corrected gradient, the watershed algorithm is used to obtain the region boundary with the extreme point obtained in S2 as the local minimum of the gradient. S4: Based on the segmented label matrix, obtain the boundary of the extreme point, perform corresponding statistical processing on the data of the region to which each extreme point belongs, extract the pit depth and diameter corresponding to the extreme point region, and statistically analyze the overall pit data distribution.

2. The method for identifying erosion pits and extracting parameters based on three-dimensional extreme points according to claim 1, characterized in that, In step S1, the original point cloud data of the corrosion specimen surface is obtained by three-dimensional laser scanning; a tilted plane is fitted by solving the plane equation of the original data points; a rotation matrix is ​​constructed based on the normal vector of the fitted plane to correct the tilt for the first time; the principal direction of the plane is determined by principal component analysis to construct a rotation matrix for the second time to align the plane with the x and y coordinate axes; and regular grid data with edges parallel to the x and y coordinate axes is generated by bilinear interpolation based on the rotated discrete plane coordinate data, and the missing values ​​of the original data are filled by nearest neighbor, finally obtaining regular and standard grid data. The specific operation steps are as follows: (1) Acquisition of raw point cloud: The raw point cloud data of the corrosion specimen surface was acquired by non-contact three-dimensional laser scanning. The original plane is tilted. (2) Sloping plane equation fitting: Assume that the plane equation constructed from the original data points is The original coordinate data of the corroded surface were obtained through non-contact three-dimensional laser scanning. Now let the original surface matrix The plane coefficients can be solved using the least squares method: ； in, This is the original height vector; (3) Rotation of the tilted plane: based on the original plane equation normal vector Normal vectors of the x, y plane Calculate the rotation axis of the two normal vectors. Rotation angle via rotating axis With rotation angle The rotation matrix can then be obtained. This corrects the plane tilt; ， ； ， ； ， in , ; Coordinates of the point after rotation ,in To obtain the original data points from the scan, the plane is then translated in the Z direction. Make the lowest point of the plane Complete the initial tilted raw data transformation to the horizontal xy plane coordinate system; (4) Second planar rotation: Extract the horizontal plane point cloud data after the first rotation and translation into two-dimensional coordinates. ; Calculate the covariance matrix of two-dimensional coordinates The eigenvalues ​​are obtained by performing eigenvalue decomposition on the covariance matrix. and the corresponding feature vectors , A rotation matrix can then be constructed. ; ； ， ； ； in , Let x and y be the mean values ​​of the point cloud in the plane, and n be the number of points. The main direction represents the direction in which the data extends the longest. This indicates that this direction is perpendicular to the main direction; Two-dimensional coordinates in the xy plane The coordinates are obtained after rotation. Align the center of the planar point cloud with the origin, align the main direction of the point cloud with the x-axis, and align the secondary direction of the point cloud with the y-axis to achieve a regularized distribution of data in the plane; (5) Grid data fitting and missing value imputation: The number of grids in the x and y directions is determined based on the x and y axis ranges of the point cloud coordinates after the second rotation and the grid spacing. ， ； ， ； ， ； ， ； in , For actively input grid precision; , The number of grid nodes in the x and y directions; Because point cloud data is discrete, each grid node height value The height will be calculated using bilinear interpolation, based on the weighted average of the four neighboring point clouds surrounding the node, where the top left point is (x1, y1, z1), the top right point is (x3, y1, z3), the bottom left point is (x1, y3, z3), and the bottom right point is (x3, y3, z4). x-direction interpolation (with fixed y): ； ； y-direction interpolation (fixed x): ； Because the initial scan data may not be rectangular or may have irregular boundaries, the obtained initial mesh boundary may not be complete. Missing values ​​will exist in the rectangular mesh matrix; these missing values ​​will be denoted as... And fill according to the nearest neighbor method, specifically: first, fill missing values ​​in the grid data along the column direction: for each column traverse row indexes ,like Then search for the distance in that column. The most recent non-missing value ,Will Fill as Then, perform a second filling along the row direction as described above; finally, obtain the vertex height data of the standard curved surface mesh. .

3. The method for identifying erosion pits and extracting parameters based on three-dimensional extreme points according to claim 1, characterized in that, In step S2, the local neighborhood range is defined by parsing control parameters. Candidate extreme points are selected based on neighborhood extreme values, mean differences, and significance. After deduplication by spatial distance constraints, the final set of maximum point coordinates of the mesh surface data is generated and used for subsequent feature analysis. The specific steps are as follows: (1) Define the control parameter: minimum significance Control the difference threshold between extreme points and their neighborhoods, minimum spacing Minimum spatial distance between extreme points, neighborhood size Define the size of a local range of rows and columns; (2) Define the local neighborhood: Let the three-dimensional surface data surface_data be a matrix Treat each testing site Calculate the neighborhood radius: , Define the neighborhood range: ； ； Where R represents a real number, and M and N are the number of rows and columns corresponding to the grid data; (3) Screening for maximum points: Denote the neighborhood submatrix Calculate the neighborhood mean for each non-boundary point. Significance : ； ； If the following conditions are met: ； Then mark For candidate maxima, their corresponding coordinates are stored in a set. Store the corresponding height value ; (4) Deduplication of Maximum Point Constraints: Spatial deduplication is performed on candidate extreme point values. according to Sort in descending order, for the sorted points Calculate and retain points Euclidean distance: ； Among them, the points have been reserved. Initially an empty set, as candidate points are calculated one by one using Euclidean distance, if... Then join in Otherwise, discard them; after deduplication, the final set of maximum points is output. .

4. The method for identifying erosion pits and extracting parameters based on three-dimensional extreme points according to claim 1, characterized in that, In step S3, after smoothing the grid surface data using Gaussian filtering, the gradient magnitude of the grid data is calculated based on the Sobel operator. The gradient magnitude is then expanded using distance transformation to determine the background and foreground regions of the gradient map. Extreme points are forcibly marked as local minima of the gradient map. Based on the modified gradient map, the watershed algorithm is applied to segment the regions and generate a label matrix, thus completing the image segmentation of the background and foreground regions. The specific steps are as follows: (1) Data smoothing: To reduce noise interference in gradient calculation, a two-dimensional Gaussian filter is used to smooth the original surface data. The smoothed data is the convolution of the original surface data and the Gaussian kernel. ; The Gaussian kernel function is: ,in, , This is the offset of a pixel within the kernel relative to the center. (2) Surface gradient calculation: The smoothed surface data is calculated using the Sobel operator. The gradient magnitude can be obtained by convolution of each gradient component: ， ； The Sobel operator kernel is defined as follows: , , The gradient magnitude can be determined based on the gradient components. calculate: ； (3) Background calculation: First, mark the spatial location of the extreme points and define the marking moments. Initialize as an all-zero matrix; for each extreme point ,set up The "background region" and the "foreground region" containing the extreme points are determined through distance transformation; the Euclidean distance from each pixel to the nearest extreme point is calculated. ； in, Let the coordinates be the k-th extreme point. When the pixel distance is greater than the threshold multiplied by the maximum distance, that is... Mark as background Update the marker simultaneously , where yuzhi is the input parameter; (4) Gradient magnitude correction: In order to constrain the "watering start point" and "region expansion range" of the watershed algorithm, and to ensure that the watershed algorithm "floods" the region from the extreme point, while avoiding the over-segmentation problem of the unconstrained watershed algorithm, the foreground extreme point and the background region will be forced to be set as local minima of the gradient topographic map through mathematical constraints, so that both become watering points at the same time. Treating the gradient magnitude as terrain height, we have: ； From the tag matrix Extract the coordinate set of all foreground extreme points: ； Where K is the total number of foreground extreme points; For each extreme point Define its domain as and all pixels in the neighborhood Perform gradient correction: ； in It is a very small positive number to ensure that the gradient value of the extreme point itself is strictly less than that of other points in the neighborhood; This represents all pixels in the neighborhood except for the extreme points; Similarly, the background region is corrected by the gradient map and forced into a local minimum region of gradient terrain, serving as the "background water injection starting point" of the watershed algorithm; From the label matrix, we can see that This is the background region, denoted as B. A region-wide gradient reduction strategy is applied to background region B to correct the gradient: ； in This is the original gradient map. Represents the minimum gradient value in the entire image. It is a very small positive number; The corrected background region becomes the flattest region globally, which can be regarded as the sea level in the watershed algorithm; The above steps complete the correction of the gradient assignment of the original data. (5) Delineation of watershed boundaries: Corrected gradient map The initial set of local minima is the union of the foreground and background regions: ； in ; The initial region and label matrix are now initialized for the subsequent watershed algorithm: Foreground area: for each initial region ; Background area: Initial area ; Tag matrix: ; in Ensure that the background label and the foreground label are not duplicated; Flood fill the corrected gradient graph based on a priority queue:

1. Water level parameters and event queue: Water level range: ,in , ; Event queue: A queue Q is used to store pixels to be processed, arranged in ascending order of gradient value h(x,y). The simulation simulates water level rising from low to high. Initially, all minimum points are added to the queue. ； 2. Expansion of the flood area: Pixels are processed in ascending order of gradient values ​​(simulating water levels rising from low to high). For each pixel (x, y), a priority queue is used. Extract and process: First, extract the pixel (x,y) with the smallest h(x,y) from Q, and let the current water level t = h(x,y); then, start neighborhood expansion based on the minimum water level, and expand the eight neighborhoods of (x,y). Pixels Perform marking and determination: If this is not marked Then It is assigned to the current region and added to the queue to await processing of its neighbors: ； If it has been marked and the tags are different Then it is marked as a pixel boundary: ； When priority queue Q is empty, all pixels are processed. At this time: Foreground area: ; Background area: ; Tag matrix: ; in This represents the foreground region excluding the boundary points. This represents the background area excluding the boundary points.

5. The method for identifying erosion pits and extracting parameters based on three-dimensional extreme points according to claim 1, characterized in that, In step S4, based on the label matrix obtained by the watershed algorithm, the set of boundary coordinates corresponding to each extreme point is first obtained; then, the depth and diameter of the pits are extracted from each extreme point region through statistical processing, thereby obtaining the pit depth features of the overall corrosion morphology. The specific steps are as follows: (1) Boundary extraction of the extreme point label matrix: Based on the label matrix obtained in S3, the boundary pixel coordinates of the label region corresponding to each extreme point are extracted and combined with the height value of the original grid data to generate three-dimensional boundary coordinates. Areas (labeled as) Extract its contour pixel coordinates. : ； (2) Pits feature extraction: Based on the contour pixel coordinates combined with the original surface height data This will generate the coordinates of the erosion pit boundary: ； Based on extreme points and boundary coordinates The depth of the erosion pit can then be determined. Diameter of the erosion pit Let the maximum and minimum values ​​of the pit boundary in the x and y directions be denoted as x, y, y, respectively. kmin x kmax y kmin y kmax So there are ; The corresponding pit depth-to-diameter ratio can be expressed as: ; (3) Distribution statistics of pits: Perform interval distribution probability statistics on the pit depth and depth-to-diameter ratio obtained in (2) to obtain the probability distribution density histogram; then fit the data to a normal distribution to obtain the corresponding data distribution.

Images