Tooth surface phase unwrapping method based on adaptive unscented information filtering

Abstract

The invention discloses a self-adaptive phase unwrapping method based on unscented information filtering. An interference pattern and a non-interference pattern are obtained through a gear tooth surface morphology error measurement system; processing the interferogram to obtain a wrapped phase diagram, and adaptively adjusting a local phase gradient window by taking the fringe density of the wrapped phase diagram as prior information; setting a threshold value according to a gray level distribution histogram of the non-interference pattern, and performing reliability judgment on pixel points in the wrapped phase diagram according to the threshold value; and unscented information filtering processing is carried out on the to-be-unwrapped pixel points, and meanwhile, the measurement noise is dynamically corrected according to the calculation result of the information, so that phase unwrapping is completed. The invention discloses a self-adaptive phase unwrapping method based on unscented information filtering, and solves the problems that an existing tooth surface phase unwrapping method is difficult to adapt to local stripe abrupt change, error accumulation of an unreliable region and low robustness.

Description

Technical Field

[0001] This invention belongs to the field of optical measurement technology, specifically relating to a tooth surface phase unwrapping method based on adaptive unscented information filtering. Background Technology

[0002] Interferometry, with its advantages of being non-contact and highly accurate, has become one of the core technologies in fields such as optical measurement, synthetic aperture radar, and medical imaging. Since the phase information carried in the interferogram needs to be obtained through phase unwrapping, its accuracy directly depends on the effectiveness of the phase unwrapping algorithm.

[0003] Currently, there are four commonly used methods for unwrapping the phase of gear tooth surfaces: path tracing, which unfolds the unwrapping path point by point, but this method is highly dependent on paths; minimum norm, which achieves phase smoothing through global optimization, but this method is sensitive to boundary conditions and requires high computational resources; deep learning, which recovers the true phase from the wrapped model by training a data model, but this method is complex and the model has a generalization problem; and interferometric image preprocessing, which preprocesses the wrapped image through filtering before applying traditional methods. All of these methods require pre-filtering the wrapped image before unwrapping, and all face the problem of over-filtering leading to the loss of some phase information or under-filtering leading to residual noise.

[0004] Therefore, phase unwrapping methods based on Bayesian filtering estimation have emerged, with Kalman filtering being one such method. Its strong noise immunity and good real-time processing capabilities make it suitable for phase unwrapping. However, the periodic structure of the involute tooth surface of gears causes the interferogram to exhibit a unique non-uniform anisotropic fringe distribution. Furthermore, the coupling noise generated by tooth surface reflection is fundamentally different from the multiplicative noise in InSAR (Interferometric Synthetic Aperture Radar), leading to limitations in the application of the aforementioned methods. Existing unwrapping methods still suffer from the following problems for gear tooth surfaces: difficulty in adapting to sudden changes in local fringe density, error accumulation in unreliable regions, and reduced robustness due to static assumptions. Summary of the Invention

[0005] The purpose of this invention is to provide a tooth surface phase unwrapping method based on adaptive unscented information filtering, which solves the problems of existing tooth surface phase unwrapping methods, such as difficulty in adapting to local stripe abrupt changes, error accumulation in unreliable regions, and low robustness.

[0006] The technical solution adopted in this invention is a tooth surface phase unwrapping method based on adaptive unscented information filtering. This method uses an interferogram and a non-interferogram obtained from a gear tooth surface topography error measurement system. The interferogram is processed to obtain a wrapped phase map. The fringe density of the wrapped phase map is used as prior information to adaptively adjust the local phase gradient window. A threshold is set based on the grayscale distribution histogram of the non-interferogram, and the reliability of pixels in the wrapped phase map is judged based on the threshold. Unscented information filtering is applied to the pixels to be unwrapped, and measurement noise is dynamically corrected based on the innovation calculation results to complete the phase unwrapping.

[0007] The invention is further characterized by: The tooth surface phase unwrapping method based on adaptive unscented information filtering specifically includes the following steps: Step 1: Using the interferogram and non-interferogram obtained by the gear tooth surface shape error measurement system, the interferogram is processed by four-step term shifting and least squares method to obtain the wrapped phase map; Step 2: Calculate the fringe density of the wrapping phase map using the cumulative gray-level difference method, standardize the fringe density, and divide the local phase gradient window according to the standardized fringe density; Step 3: Adaptively calculate the local phase gradient estimate and the local phase gradient estimation error based on the local phase gradient window; Step 4: Calculate the quality chart; Step 5: Set a threshold based on the gray-level distribution histogram of the non-interferometric image, and divide the wrapped phase image into reliable and unreliable regions based on the threshold; Step 6: Select the point with the best quality in the quality map within the reliable region as the seed point. Starting from the seed point, select the point with the best quality in the neighborhood of the seed point and within the reliable region as the pixel to be unwrapped. Step 7: Perform unscented information filtering to solve the information matrix and information vector of the unwrapped pixel, and solve the measurement vector covariance matrix and cross-covariance matrix of the measurement vector of the unwrapped pixel. Step 8: Calculate the new information and adaptively adjust the measurement noise based on the new information; Step 9: Combine the results of Step 7 and Step 8 to perform phase unwrapping on the pixels to be wrapped; Step 10: Select the highest quality pixel within the neighborhood of the unwrapped pixel and within the reliable region as the next pixel to be processed; repeat steps 7-9 until all pixels in the reliable region have been unwrapped; Step 11: Select the point with the best quality in the quality map in the unreliable region as the seed point. Starting from the seed point, select the point with the best quality in the neighborhood of the seed point that is in the unreliable region as the pixel to be unwrapped. Unwrap the pixel according to steps 7 to 9. Select the point with the best quality in the neighborhood of the unwrapped pixel that is in the unreliable region as the next point to be processed. Repeat steps 7 to 9 until all pixels in the unreliable region are unwrapped.

[0008] In step 2, the fringe density of the wrapped phase map is calculated using the cumulative gray-level difference method. The calculation formula is as follows: (1) In the formula: , These represent the cumulative grayscale difference along the x and y directions, respectively. Represents pixels The actual grayscale value is shown; B represents the stripe width; L represents the window size. This represents the desired fringe density; The stripe density is standardized and converted into a normal distribution with a mean of 0 and a standard deviation of 1, as shown in the following formula: (2) In the formula, Indicates stripe density, Indicates the mean value of the stripes. The standard deviation of the stripe density; Based on the standardized fringe density, the local phase gradient estimation window is... It is divided into three parts, as shown below: (3).

[0009] Step 3 specifically includes the following sub-steps: Step 3.1: Substitute the wrapped phase into the exponential function Among them To enclose the phase, a complex interferogram containing phase information is obtained; Step 3.2: Select pixels in the complex interferogram Centered on, the window is The matrix is ​​set as ; Step 3.3: Perform singular value decomposition on the expression, as shown below: (4) In the formula, yes The calculation results obtained after singular value decomposition are all Matrix; in, (5) In the formula, Principal singular values, It is a non-principal singular value; Formula (4) can be rewritten as follows: (6); Step 3.4: Extract the three matrices , means as follows: (7) Step 3.5: For the matrix Perform singular value decomposition to obtain principal singular values. Left principal singular vector and right principal singular vector , means as follows: (8) In the formula, , , Indicates the remaining components besides the main components; Formula (7) can be rewritten as follows: (9) Pixels The local phase gradients in the row and column directions are obtained using the following formulas: (10) In the formula, , Local phase gradients in the row and column directions, They are respectively The pseudo-inverse matrix, Indicates the calculation of the argument. This indicates how to find the conjugate of a complex number; In a local window, pixels relative adjacent pixels The phase gradient estimate in the direction is calculated as follows: (11) In the formula, Represents pixels relative adjacent pixels Phase gradient in the direction, , Representing adjacent pixels and pixels The phase value; Step 3.6: Calculate the local phase gradient estimation error, as shown below: (12) In the formula, , These represent the lengths of the local window in the y and x directions, respectively. For pixels The pseudo-coherence coefficient; This indicates the calculation of the error between two adjacent pixels; Represents pixels relative adjacent pixels The local phase gradient estimation error.

[0010] Step 4 specifically involves selecting any one of the following as a parameter to measure the quality of the wrapped phase map: pseudo-coherence coefficient, phase derivative variance, or maximum phase gradient. The quality of each pixel is then calculated using the parameter formula to construct the quality map. The formula for calculating the pseudocoherence coefficient is as follows: (13) In the formula, Represents pixels The pseudocoherence coefficient, Represents pixels The wrapping phase value, k represents the size of the neighborhood window centered on the target pixel; The formula for calculating the variance of the phase derivative is as follows: (14) In the formula, Represents pixels The phase derivative variance This represents the partial derivative of the wrapping phase in the x-direction; This represents the partial derivative of the wrapping phase in the y-direction; express exist The average value within the window; express exist The average value within the window; where, , , This indicates that the calculation result is restricted to a specific package operator; The formula for calculating the maximum phase gradient is as follows: (15) In the formula, Represents pixels The maximum phase gradient.

[0011] Step 5, the method for setting the threshold, specifically includes the following sub-steps: Step 5.1: Calculate the residual points in the package phase, construct a gray value histogram, and record the location of the residual points and the corresponding gray values; Step 5.2: Count the number of residual points corresponding to each gray value, and establish the relationship between the percentage of residual points and the gray value; Step 5.3: From the percentage of residual points Starting from 0, according to The amplitude increases incrementally, and the grayscale value is recorded. ; Step 5.4: Stop increasing the threshold when the following formula condition is met. The gray value corresponding to the residual percentage at this point is the threshold: (16) In the formula, k Indicates the number of times; The non-interferometric image is divided according to the following conditions: the region composed of pixels with gray values ​​greater than or equal to the threshold is the reliable region, and the region composed of pixels with gray values ​​less than the threshold is the unreliable region; at the same time, since the non-interferometric image corresponds one-to-one with the pixels in the wrapped phase image, the wrapped phase image is also divided into reliable and unreliable regions.

[0012] Step 7 specifically includes: The pixels to be unpacked The initial values ​​of the state variables and state covariance matrix are estimated by incorporating them into the unscented information filter, as shown below: (17) In the formula, Represents state variables The initial value; Represents the state covariance matrix The initial value; Indicates pixels that have been unwrapped. The weights; This indicates that the unwrapped pixels have been identified. State variables; This indicates that the unwrapped pixels have been identified. The covariance matrix; Represents pixels The corresponding covariance matrix; Represents pixels Relative to pixels The local phase gradient estimate; Indicates the signal-to-noise ratio; Represents pixels Signal-to-noise ratio; Represents pixels The set of unwrapped pixels within the neighborhood; Initialize the state variables Initial values ​​of the state covariance matrix The corresponding Sigma point set is obtained by performing an unscented transformation, as follows: (18) In the formula, Represents the Sigma point. d This represents the dimension of the state variables. The scaling parameter of the Sigma point is represented by, where and For tuning factor, Take 0, ; The weights corresponding to each Sigma point are calculated as follows: (19) In the formula, They represent The weights of the state variables and the state covariance matrix; Indicates respectively The weights of the state variables and the state covariance matrix; Indicates the regulating factor. Take 2; Based on the calculated Sigma point set, for each pixel... The prediction and estimation of the state variables, as well as the corresponding prediction and estimation of the state covariance matrix, are calculated as follows: (20) In the formula, This represents the predicted estimate of the state variable. This represents the predicted estimate of the state covariance matrix; Based on the definition of information filtering, calculate the information matrix. and information vector , means as follows: (twenty one) For pixels To predict the measurement vector, substitute the Sigma point set into the measurement equation and calculate the measurement vector covariance matrix and cross-covariance matrix. The calculation formulas are as follows: (twenty two) In the formula, This represents the measurement value of each point calculated using the Sigma point set; This indicates that pixel points are estimated using the Sigma point set. The measured value, i.e. the predicted measured value; For the measurement equation, ; The cross-covariance matrix representing the measurement vectors; This represents the covariance matrix of the measurement vector.

[0013] Step 8 specifically includes: The innovation is the difference between the actual measurement and the predicted measurement, expressed as follows: (twenty three) In the formula, Indicates new information; Indicates the measured value; Indicates the predicted measurement value; Based on the Mahalanobis distance and the statistics of the information obtained during the filtering process, and given that they follow a chi-square distribution, it is expressed as follows: (twenty four) In the formula, Indicates the transposition of new information; Indicates the measurement noise covariance; A statistical measure representing new information; For the dimension of the new information, This indicates that the distribution follows a chi-square distribution; The threshold is determined by the chi-square distribution of the new information, and the measurement noise adjustment coefficient is calculated as follows: (25) (26) In the formula, Indicates the measurement noise adjustment factor; Represents pixels The amount of new information statistics; This represents the updated measurement noise covariance; , , Represents pixels Signal-to-noise ratio at the location.

[0014] Step 9 specifically includes: The information contribution is calculated using the following formula: (27) In the formula, Observational contribution items representing information state Represents the observation contribution term of the information matrix; For pixels Phase expansion is performed, as follows: (28) In the formula, , , Representing pixels The updated information state, the updated information matrix, and the final covariance matrix. Represents pixels The unfolded phase.

[0015] The beneficial effects of this invention are: This invention presents a tooth surface phase unwrapping method based on adaptive unscented information filtering (UIF). It adaptively adjusts the window size in the phase gradient using the fringe density of the wrapped phase map as prior information, solving the problem that a fixed window cannot accurately estimate the magnitude of the local phase gradient and improving the accuracy of phase estimation. It distinguishes reliable and unreliable regions based on grayscale values, optimizes the unwrapping path, and reduces error accumulation in unreliable regions, resulting in more accurate unwrapping results. Furthermore, it dynamically corrects measurement noise based on the innovation calculation results, reducing the root mean square error and making the unwrapping results more continuous, smooth, and stable. Compared to traditional methods and UIF, the method used to process gear tooth surface wrapped phase maps reduces the number of discontinuous phase points by 25-35%. Attached Figure Description

[0016] Figure 1 This is the wrapped phase diagram from the simulation experiment of the tooth surface phase unwrapping method based on adaptive unscented information filtering of this invention; Figure 2 The fringe density map in the simulation experiment of the tooth surface phase unwrapping method based on adaptive unscented information filtering in this invention; Figure 3 The simulation results show a comparison of the unwrapping results and discontinuous phase point results of this invention with other methods. Figure 4 This is a comparison chart of the unwrapping curves of the present invention and other methods in the simulation experiment for the 110th row. Detailed Implementation

[0017] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments.

[0018] This invention relates to a tooth surface phase unwrapping method based on adaptive unscented information filtering, which specifically includes the following steps: Step 1: Using the interferogram and non-interferogram obtained by the gear tooth surface shape error measurement system, the interferogram is processed by four-step term shifting and least squares method to obtain the wrapped phase map.

[0019] The four-step phase-shifting method is a commonly used technique for extracting the wrapped phase from interferograms. Its core is to acquire four interferograms with fixed phase differences (typically 0, π / 2, π, and 3π / 2) and use this phase difference information to calculate the original wrapped phase. After obtaining the wrapped phase, the least squares method is often used for phase gradient estimation or initial phase unfolding, laying the foundation for subsequent accurate unwrapping.

[0020] Step 2: Calculate the fringe density of the wrapping phase map using the cumulative gray-level difference method, standardize the fringe density, and divide the local phase gradient window according to the standardized fringe density.

[0021] The formula for calculating stripe density is as follows: (1) In the formula: , These represent the cumulative grayscale difference along the x and y directions, respectively. Represents pixels The actual grayscale value is shown; B represents the stripe width; L represents the window size. This represents the desired fringe density; The stripe density is standardized and converted into a normal distribution with a mean of 0 and a standard deviation of 1, as shown in the following formula: (2) In the formula, Indicates stripe density, Indicates the mean value of the stripes. The standard deviation of the stripe density.

[0022] Because of the non-uniform distribution of gear stripes, a fixed window cannot accurately estimate the magnitude of the local phase gradient in local phase gradient estimation. Dynamically adjusting the window size based on the standardized stripe density—smaller window for higher density, larger window for lower density—leads to more accurate calculations. The local phase gradient estimation window is adjusted according to the standardized stripe density. It is divided into three parts, as shown below: (3).

[0023] Step 3: Adaptively calculate the local phase gradient estimate and the local phase gradient estimation error based on the local phase gradient window. This includes the following sub-steps: Step 3.1: Substitute the wrapped phase into the exponential function Among them To enclose the phase, a complex interferogram containing phase information is obtained; Step 3.2: Select pixels in the complex interferogram Centered on, the window is The matrix is ​​set as ; Step 3.3: Perform singular value decomposition on the expression, as shown below: (4) In the formula, yes The calculation results obtained after singular value decomposition are all Matrix; in, (5) In the formula, Principal singular values, It is a non-principal singular value; Formula (4) can be rewritten as follows: (6); Step 3.4: Extract the three matrices , means as follows: (7) Step 3.5: For the matrix Perform singular value decomposition to obtain principal singular values. Left principal singular vector and right principal singular vector , means as follows: (8) In the formula, , , Indicates the remaining components besides the main components; Formula (7) can be rewritten as follows: (9) Pixels The local phase gradients in the row and column directions are obtained using the following formulas: (10) In the formula, , Local phase gradients in the row and column directions, They are respectively The pseudo-inverse matrix, Indicates the calculation of the argument. This indicates how to find the conjugate of a complex number; In a local window, pixels relative adjacent pixels The phase gradient estimate in the direction is calculated as follows: (11) In the formula, Represents pixels relative adjacent pixels Phase gradient in the direction, , Representing adjacent pixels and pixels The phase value; Step 3.6: Calculate the local phase gradient estimation error, as shown below: (12) In the formula, , These represent the lengths of the local window in the y and x directions, respectively. For pixels The pseudo-coherence coefficient; This indicates the calculation of the error between two adjacent pixels; Represents pixels relative adjacent pixels The local phase gradient estimation error.

[0024] Step 4: Calculate the quality chart; A quality map is used to measure the quality of each pixel in a wrapper phase map. Commonly used parameters to reflect phase quality include pseudocoherence diff (PSD), phase derivative deviation (PDV), and maximum phase gradient (MPG). The formulas for calculating these parameters are as follows: The formula for calculating the pseudocoherence coefficient is as follows: (13) In the formula, Represents pixels The pseudocoherence coefficient, Represents pixels The wrapping phase value, k represents the size of the neighborhood window centered on the target pixel; The formula for calculating the variance of the phase derivative is as follows: (14) In the formula, Represents pixels The phase derivative variance This represents the partial derivative of the wrapping phase in the x-direction; This represents the partial derivative of the wrapping phase in the y-direction; express exist The average value within the window; express exist The average value within the window; where, , , This indicates that the calculation result is restricted to a specific package operator; The formula for calculating the maximum phase gradient is as follows: (15) In the formula, Represents pixels The maximum phase gradient.

[0025] Choose one of the three parameters mentioned above as the parameter for measuring the quality of the phase map, and calculate the quality of each pixel according to the corresponding formula to construct a quality map. The quality map provides path support for the following pixel unwrapping process.

[0026] Step 5: Set a threshold based on the grayscale distribution histogram of the non-interferometric image, and divide the wrapped phase image into reliable and unreliable regions based on the threshold.

[0027] The grayscale distribution in the non-interference image of the tooth surface exhibits abrupt changes due to variations in light intensity distribution on the CCD caused by the actual use and processing of the gear. Regions with lower grayscale values ​​in the non-interference image of the tooth surface are caused by weaker signal light intensity, making them susceptible to noise during phase unwrapping and leading to unreliable results. To prevent repeated threshold adjustments, an adaptive threshold adjustment method is used to set the threshold based on the grayscale distribution histogram. The specific steps for threshold setting are as follows: Step 5.1: Calculate the residual points in the package phase, construct a gray value histogram, and record the location of the residual points and the corresponding gray values; Step 5.2: Count the number of residual points corresponding to each gray value, and establish the relationship between the percentage of residual points and the gray value; Step 5.3: From the percentage of residual points Starting from 0, according to The amplitude increases incrementally, and the grayscale value is recorded. ; Step 5.4: Stop increasing the threshold when the following formula condition is met. The gray value corresponding to the residual percentage at this point is the threshold: (16) In the formula, k Indicates the number of times; After calculating the threshold, the non-interference map is divided according to the following conditions: the region composed of pixels with gray values ​​greater than or equal to the threshold is the reliable region, and the region composed of pixels with gray values ​​less than the threshold is the unreliable region; at the same time, since the non-interference map of the tooth surface corresponds one-to-one with the pixels in the wrapped phase map, the wrapped phase map is also divided into reliable and unreliable regions.

[0028] Step 6: Select the point with the best quality in the quality map within the reliable region as the seed point. Starting from the seed point, select the point with the best quality in the neighborhood of the seed point and within the reliable region as the pixel to be unwrapped.

[0029] Step 7: Perform unscented information filtering to solve for the information matrix and information vector of the pixels to be unwrapped, and solve for the covariance matrix and cross-covariance matrix of the measurement vectors of the pixels to be unwrapped. The specific process is as follows: The pixels to be unpacked The initial values ​​of the state variables and state covariance matrix are estimated by incorporating them into the unscented information filter, as shown below: (17) In the formula, Represents state variables The initial value; Represents the state covariance matrix The initial value; Indicates pixels that have been unwrapped. The weights; This indicates that the unwrapped pixels have been identified. State variables; This indicates that the unwrapped pixels have been identified. The covariance matrix; Represents pixels The corresponding covariance matrix; Represents pixels Relative to pixels The local phase gradient estimate; Indicates the signal-to-noise ratio; Represents pixels Signal-to-noise ratio; Represents pixels The set of unwrapped pixels within the neighborhood; Initialize the state variables Initial values ​​of the state covariance matrix The corresponding Sigma point set is obtained by performing an unscented transformation, as follows: (18) In the formula, Represents the Sigma point. d This represents the dimension of the state variables. The scaling parameter of the Sigma point is represented by, where and For tuning factor, Take 0, ; The weights corresponding to each Sigma point are calculated as follows: (19) In the formula, They represent The weights of the state variables and the state covariance matrix; Indicates respectively The weights of the state variables and the state covariance matrix; Indicates the regulating factor. Take 2; Based on the calculated Sigma point set, for each pixel... The prediction and estimation of the state variables, as well as the corresponding prediction and estimation of the state covariance matrix, are calculated as follows: (20) In the formula, This represents the predicted estimate of the state variable. This represents the predicted estimate of the state covariance matrix; Based on the definition of information filtering, calculate the information matrix. and information vector , means as follows: (twenty one) For pixels To predict the measurement vector, substitute the Sigma point set into the measurement equation and calculate the measurement vector covariance matrix and cross-covariance matrix. The calculation formulas are as follows: (twenty two) In the formula, This represents the measurement value of each point calculated using the Sigma point set; This indicates that pixel points are estimated using the Sigma point set. The measured value, i.e. the predicted measured value; For the measurement equation, ; The cross-covariance matrix representing the measurement vectors; This represents the covariance matrix of the measurement vector.

[0030] Step 8: Calculate the new information and adaptively adjust the measurement noise based on the new information.

[0031] The innovation is the difference between the actual measurement and the predicted measurement, expressed as follows: (twenty three) In the formula, Indicates new information; Indicates the measured value; Indicates the predicted measurement value; Based on the Mahalanobis distance and the statistics of the information obtained during the filtering process, and given that they follow a chi-square distribution, it is expressed as follows: (twenty four) In the formula, Indicates the transposition of new information; Indicates the measurement noise covariance; A statistical measure representing new information; For the dimension of the new information, This indicates that the distribution follows a chi-square distribution; Mahalanobis distance, also known as Mahalanobis distance, is a distance metric used to measure the relationship between a data point and the center of a distribution, and can measure the degree of data anomalies. The newly developed Mahalanobis distance measures the normalization of observation errors and can measure measurement noise; excessive or insufficient measurement noise will affect filtering performance. The accuracy of the unpacking results is achieved by dynamically adjusting the measurement noise.

[0032] The threshold is determined by the chi-square distribution of the new information, and the measurement noise adjustment coefficient is calculated as follows: (25) (26) In the formula, Indicates the measurement noise adjustment factor; Represents pixels The amount of new information statistics; This represents the updated measurement noise covariance; , , Represents pixels Signal-to-noise ratio at the location.

[0033] Step 9: Combine the results of Step 7 and Step 8 to perform phase unwrapping on the pixels to be wrapped.

[0034] The information contribution is calculated using the following formula: (27) In the formula, Observational contribution items representing information state Represents the observation contribution term of the information matrix; For pixels Phase expansion is performed, as follows: (28) In the formula, , , Representing pixels The updated information state, the updated information matrix, and the final covariance matrix. Represents pixels The unfolded phase.

[0035] Step 10: Select the highest quality pixel within the neighborhood of the unwrapped pixel and within the reliable region as the next pixel to be processed; repeat steps 7-9 until all pixels in the reliable region have been unwrapped; Step 11: Select the point with the best quality in the quality map in the unreliable region as the seed point. Starting from the seed point, select the point with the best quality in the neighborhood of the seed point that is in the unreliable region as the pixel to be unwrapped. Unwrap the pixel according to steps 7 to 9. Select the point with the best quality in the neighborhood of the unwrapped pixel that is in the unreliable region as the next point to be processed. Repeat steps 7 to 9 until all pixels in the unreliable region are unwrapped.

[0036] The method of this invention combines unscented information filtering phase unwrapping based on local phase gradient estimation with stripe density features, and realizes a gear tooth surface phase unwrapping method based on adaptive unscented information filtering based on phase reliability judgment of non-interferometric images of the tooth surface and adaptive measurement noise.

[0037] Example 1 This embodiment provides a tooth surface phase unwrapping method based on adaptive unscented information filtering, which specifically includes the following steps: Step 1: Using the interferogram and non-interferogram obtained by the gear tooth surface shape error measurement system, the interferogram is processed by four-step term shifting and least squares method to obtain the wrapped phase map; Step 2: Calculate the fringe density of the wrapping phase map using the cumulative gray-level difference method, standardize the fringe density, and divide the local phase gradient window according to the standardized fringe density; Step 3: Adaptively calculate the local phase gradient estimate and the local phase gradient estimation error based on the local phase gradient window; Step 4: Calculate the quality chart; Step 5: Set a threshold based on the gray-level distribution histogram of the non-interferometric image, and divide the wrapped phase image into reliable and unreliable regions based on the threshold; Step 6: Select the point with the best quality in the quality map within the reliable region as the seed point. Starting from the seed point, select the point with the best quality in the neighborhood of the seed point and within the reliable region as the pixel to be unwrapped. Step 7: Perform unscented information filtering to solve the information matrix and information vector of the unwrapped pixel, and solve the measurement vector covariance matrix and cross-covariance matrix of the measurement vector of the unwrapped pixel. Step 8: Calculate the new information and adaptively adjust the measurement noise based on the new information; Step 9: Combine the results of Step 7 and Step 8 to perform phase unwrapping on the pixels to be wrapped; Step 10: Select the highest quality pixel within the neighborhood of the unwrapped pixel and within the reliable region as the next pixel to be processed; repeat steps 7-9 until all pixels in the reliable region have been unwrapped; Step 11: Select the point with the best quality in the quality map in the unreliable region as the seed point. Starting from the seed point, select the point with the best quality in the neighborhood of the seed point that is in the unreliable region as the pixel to be unwrapped. Unwrap the pixel according to steps 7 to 9. Select the point with the best quality in the neighborhood of the unwrapped pixel that is in the unreliable region as the next point to be processed. Repeat steps 7 to 9 until all pixels in the unreliable region are unwrapped.

[0038] Example 2 Based on Example 1, in step 2, the fringe density of the wrapped phase map is calculated using the cumulative gray-scale difference method. The calculation formula is as follows: (1) In the formula: , These represent the cumulative grayscale difference along the x and y directions, respectively. Represents pixels The actual grayscale value is shown; B represents the stripe width; L represents the window size. This represents the desired fringe density; The stripe density is standardized and converted into a normal distribution with a mean of 0 and a standard deviation of 1, as shown in the following formula: (2) In the formula, Indicates stripe density, Indicates the mean value of the stripes. The standard deviation of the stripe density; Based on the standardized fringe density, the local phase gradient estimation window is... It is divided into three parts, as shown below: (3).

[0039] Example 3 Based on Example 2, step 3 specifically includes the following sub-steps: Step 3.1: Substitute the wrapped phase into the exponential function Among them To enclose the phase, a complex interferogram containing phase information is obtained; Step 3.2: Select pixels in the complex interferogram Centered on, the window is The matrix is ​​set as ; Step 3.3: Perform singular value decomposition on the expression, as shown below: (4) In the formula, yes The calculation results obtained after singular value decomposition are all Matrix; in, (5) In the formula, Principal singular values, It is a non-principal singular value; Formula (4) can be rewritten as follows: (6); Step 3.4: Extract the three matrices , means as follows: (7) Step 3.5: For the matrix Perform singular value decomposition to obtain principal singular values. Left principal singular vector and right principal singular vector , means as follows: (8) In the formula, , , Indicates the remaining components besides the main components; Formula (7) can be rewritten as follows: (9) Pixels The local phase gradients in the row and column directions are obtained using the following formulas: (10) In the formula, , Local phase gradients in the row and column directions, They are respectively The pseudo-inverse matrix, Indicates the calculation of the argument. This indicates how to find the conjugate of a complex number; In a local window, pixels relative adjacent pixels The phase gradient estimate in the direction is calculated as follows: (11) In the formula, Represents pixels relative adjacent pixels Phase gradient in the direction, , Representing adjacent pixels and pixels The phase value; Step 3.6: Calculate the local phase gradient estimation error, as shown below: (12) In the formula, , These represent the lengths of the local window in the y and x directions, respectively. For pixels The pseudo-coherence coefficient; This indicates the calculation of the error between two adjacent pixels; Represents pixels relative adjacent pixels The local phase gradient estimation error.

[0040] Example 4 Based on Example 3, step 4 specifically involves: selecting any one of the pseudo-coherence coefficient, phase derivative variance, and maximum phase gradient as a parameter to measure the quality of the wrapped phase map, calculating the quality of each pixel according to the parameter formula, and constructing a quality map. The formula for calculating the pseudocoherence coefficient is as follows: (13) In the formula, Represents pixels The pseudocoherence coefficient, Represents pixels The wrapping phase value, k represents the size of the neighborhood window centered on the target pixel; The formula for calculating the variance of the phase derivative is as follows: (14) In the formula, Represents pixels The phase derivative variance This represents the partial derivative of the wrapping phase in the x-direction; This represents the partial derivative of the wrapping phase in the y-direction; express exist The average value within the window; express exist The average value within the window; where, , , This indicates that the calculation result is restricted to a specific package operator; The formula for calculating the maximum phase gradient is as follows: (15) In the formula, Represents pixels The maximum phase gradient.

[0041] Step 5, the method for setting the threshold, specifically includes the following sub-steps: Step 5.1: Calculate the residual points in the package phase, construct a gray value histogram, and record the location of the residual points and the corresponding gray values; Step 5.2: Count the number of residual points corresponding to each gray value, and establish the relationship between the percentage of residual points and the gray value; Step 5.3: From the percentage of residual points Starting from 0, according to The amplitude increases incrementally, and the grayscale value is recorded. ; Step 5.4: Stop increasing the threshold when the following formula condition is met. The gray value corresponding to the residual percentage at this point is the threshold: (16) In the formula, k Indicates the number of times; The non-interferometric image is divided according to the following conditions: the region composed of pixels with gray values ​​greater than or equal to the threshold is the reliable region, and the region composed of pixels with gray values ​​less than the threshold is the unreliable region; at the same time, since the non-interferometric image corresponds one-to-one with the pixels in the wrapped phase image, the wrapped phase image is also divided into reliable and unreliable regions.

[0042] Example 5 Based on Example 4, step 7 specifically includes: The pixels to be unpacked The initial values ​​of the state variables and state covariance matrix are estimated by incorporating them into the unscented information filter, as shown below: (17) In the formula, Represents state variables The initial value; Represents the state covariance matrix The initial value; Indicates pixels that have been unwrapped. The weights; This indicates that the unwrapped pixels have been identified. State variables; This indicates that the unwrapped pixels have been identified. The covariance matrix; Represents pixels The corresponding covariance matrix; Represents pixels Relative to pixels The local phase gradient estimate; Indicates the signal-to-noise ratio; Represents pixels Signal-to-noise ratio; Represents pixels The set of unwrapped pixels within the neighborhood; Initialize the state variables Initial values ​​of the state covariance matrix The corresponding Sigma point set is obtained by performing an unscented transformation, as follows: (18) In the formula, Represents the Sigma point. d This represents the dimension of the state variables. The scaling parameter of the Sigma point is represented by, where and For tuning factor, Take 0, ; The weights corresponding to each Sigma point are calculated as follows: (19) In the formula, They represent The weights of the state variables and the state covariance matrix; Indicates respectively The weights of the state variables and the state covariance matrix; Indicates the regulating factor. Take 2; Based on the calculated Sigma point set, for each pixel... The prediction and estimation of the state variables, as well as the corresponding prediction and estimation of the state covariance matrix, are calculated as follows: (20) In the formula, This represents the predicted estimate of the state variable. This represents the predicted estimate of the state covariance matrix; Based on the definition of information filtering, calculate the information matrix. and information vector , means as follows: (twenty one) For pixels To predict the measurement vector, substitute the Sigma point set into the measurement equation and calculate the measurement vector covariance matrix and cross-covariance matrix. The calculation formulas are as follows: (twenty two) In the formula, This represents the measurement value of each point calculated using the Sigma point set; This indicates that pixel points are estimated using the Sigma point set. The measured value, i.e. the predicted measured value; For the measurement equation, ; The cross-covariance matrix representing the measurement vectors; This represents the covariance matrix of the measurement vector.

[0043] Step 8 specifically includes: The innovation is the difference between the actual measurement and the predicted measurement, expressed as follows: (twenty three) In the formula, Indicates new information; Indicates the measured value; Indicates the predicted measurement value; Based on the Mahalanobis distance and the statistics of the information obtained during the filtering process, and given that they follow a chi-square distribution, it is expressed as follows: (twenty four) In the formula, Indicates the transposition of new information; Indicates the measurement noise covariance; A statistical measure representing new information; For the dimension of the new information, This indicates that the distribution follows a chi-square distribution; The threshold is determined by the chi-square distribution of the new information, and the measurement noise adjustment coefficient is calculated as follows: (25) (26) In the formula, Indicates the measurement noise adjustment factor; Represents pixels The amount of new information statistics; This represents the updated measurement noise covariance; , , Represents pixels Signal-to-noise ratio at the location.

[0044] Example 6 Based on Example 5, step 9 specifically includes: The information contribution is calculated using the following formula: (27) In the formula, Observational contribution items representing information state Represents the observation contribution term of the information matrix; For pixels Phase expansion is performed, as follows: (28) In the formula, , , Representing pixels The updated information state, the updated information matrix, and the final covariance matrix. Represents pixels The unfolded phase.

[0045] Simulation Experiment The simulation was run on a computer in the MATLAB (MATLAB R2022a) environment, equipped with an Intel(R) Core(TM) i7-8550U + 64-bit Windows 10 Professional Edition + 4GB RAM.

[0046] The wrapped phase map generated using the peaks function with a coefficient of 6 and a signal-to-noise ratio of 1.9501 dB is as follows: Figure 1 As shown, and after calculating the fringe density, it forms as shown. Figure 2 The fringe density map is shown. Traditional unwrapping methods—least squares, quality map method, branching method, and unscented information filtering method—are used to analyze it. Figure 1 The unwrapping is performed using the wrapped phase map. The unwrapping results of the method of this invention and the traditional method are shown in the figure below. Figure 3 As shown, where Figure 3 Images (a) through (e) show the unwrapping results obtained using the Least Squares Method (LSM), Quality Guide Phase Unwrapping (QGPU), Branch Cut Algorithm (BUT), Unscented Information Filtering (UIF), and the method of this invention, respectively. Figure 3 In the diagram, (f) to (g) represent the discontinuous phase point results obtained using the above method. Figure 3 As can be seen from Table 1, the method of this invention results in the fewest discontinuous phase points after unwrapping, a decrease of 47.9% compared to the UIF method. Figure 4 To use four conventional methods and the method of the present invention to Figure 1 The comparison chart of unpacking results in row 110 shows that the unpacking result obtained by the method of this invention is the smoothest.

[0047] Table 1 Simulation Experiment Data Results Parameters

Claims

1. A tooth surface phase unwrapping method based on adaptive unscented information filtering, characterized in that, Interferograms and non-interferograms are obtained through a gear tooth surface shape error measurement system; the interferogram is processed to obtain a wrapped phase map, and the local phase gradient window is adaptively adjusted using the fringe density of the wrapped phase map as prior information; a threshold is set according to the gray-level distribution histogram of the non-interferogram, and the reliability of the pixels in the wrapped phase map is judged according to the threshold; the unscented information filtering is performed on the pixels to be unwrapped, and the measurement noise is dynamically corrected according to the innovation calculation results to complete the phase unwrapping.

2. The tooth surface phase unwrapping method based on adaptive unscented information filtering according to claim 1, characterized in that, Specifically, the steps include the following: Step 1: Using the interferogram and non-interferogram obtained by the gear tooth surface shape error measurement system, the interferogram is processed by four-step term shifting and least squares method to obtain the wrapped phase map; Step 2: Calculate the fringe density of the wrapping phase map using the cumulative gray-level difference method, standardize the fringe density, and divide the local phase gradient window according to the standardized fringe density; Step 3: Adaptively calculate the local phase gradient estimate and the local phase gradient estimation error based on the local phase gradient window; Step 4: Calculate the quality chart; Step 5: Set a threshold based on the gray-level distribution histogram of the non-interferometric image, and divide the wrapped phase image into reliable and unreliable regions based on the threshold; Step 6: Select the point with the best quality in the quality map within the reliable region as the seed point. Starting from the seed point, select the point with the best quality in the neighborhood of the seed point and within the reliable region as the pixel to be unwrapped. Step 7: Perform unscented information filtering to solve the information matrix and information vector of the unwrapped pixel, and solve the measurement vector covariance matrix and cross-covariance matrix of the measurement vector of the unwrapped pixel. Step 8: Calculate the new information and adaptively adjust the measurement noise based on the new information; Step 9: Combine the results of Step 7 and Step 8 to perform phase unwrapping on the pixels to be wrapped; Step 10: Select the highest quality pixel within the neighborhood of the unwrapped pixel and within the reliable region as the next pixel to be processed; repeat steps 7-9 until all pixels in the reliable region have been unwrapped; Step 11: Select the point with the best quality in the quality map in the unreliable region as the seed point. Starting from the seed point, select the point with the best quality in the neighborhood of the seed point that is in the unreliable region as the pixel to be unwrapped. Unwrap the pixel according to steps 7 to 9. Select the point with the best quality in the neighborhood of the unwrapped pixel that is in the unreliable region as the next point to be processed. Repeat steps 7 to 9 until all pixels in the unreliable region are unwrapped.

3. The tooth surface phase unwrapping method based on adaptive unscented information filtering according to claim 2, characterized in that, In step 2, the fringe density of the wrapped phase map is calculated using the cumulative gray-level difference method. The calculation formula is as follows: (1) In the formula: , These represent the cumulative grayscale difference along the x and y directions, respectively. Represents pixels The actual grayscale value is shown; B represents the stripe width; L represents the window size. This represents the desired fringe density; The stripe density is standardized and converted into a normal distribution with a mean of 0 and a standard deviation of 1, as shown in the following formula: （2） In the formula, Indicates stripe density, Indicates the average value of the stripes. The standard deviation of the stripe density; Based on the standardized fringe density, the local phase gradient estimation window is... It is divided into three parts, as shown below: （3）。 4. The tooth surface phase unwrapping method based on adaptive unscented information filtering according to claim 2, characterized in that, Step 3 specifically includes the following sub-steps: Step 3.1: Substitute the wrapped phase into the exponential function Among them To encapsulate the phase, a complex interferogram containing phase information is obtained; Step 3.2: Select pixels in the complex interferogram Centered on, the window is The matrix is ​​set as ; Step 3.3: Perform singular value decomposition on the expression, as shown below: （4） In the formula, yes The calculation results obtained after singular value decomposition are all Matrix; in, （5） In the formula, Principal singular values, It is a non-principal singular value; Formula (4) can be rewritten as follows: （6）； Step 3.4: Extract the three matrices , means as follows: （7） Step 3.5: For the matrix Perform singular value decomposition to obtain principal singular values. Left principal singular vector and right principal singular vector , means as follows: （8） In the formula, , , Indicates the remaining components besides the main components; Formula (7) can be rewritten as follows: （9） Pixels The local phase gradients in the row and column directions are obtained using the following formulas: （10） In the formula, , Local phase gradients in the row and column directions, They are respectively The pseudo-inverse matrix, Indicates the calculation of the argument. This indicates how to find the conjugate of a complex number; In a local window, pixels relative adjacent pixels The phase gradient estimate in the direction is calculated as follows: （11） In the formula, Represents pixels relative adjacent pixels Phase gradient in the direction, , Representing adjacent pixels and pixels The phase value; Step 3.6: Calculate the local phase gradient estimation error, as shown below: （12） In the formula, , These represent the lengths of the local window in the y and x directions, respectively. For pixels The pseudo-coherence coefficient; This indicates the calculation of the error between two adjacent pixels; Represents pixels relative adjacent pixels The local phase gradient estimation error.

5. The tooth surface phase unwrapping method based on adaptive unscented information filtering according to claim 2, characterized in that, Step 4 specifically involves selecting any one of the following as a parameter to measure the quality of the wrapped phase map: pseudo-coherence coefficient, phase derivative variance, or maximum phase gradient. The quality of each pixel is then calculated using the parameter formula to construct the quality map. The formula for calculating the pseudo-coherence coefficient is as follows: (13) In the formula, Represents pixels The pseudocoherence coefficient, Represents pixels The wrapping phase value, k represents the size of the neighborhood window centered on the target pixel; The formula for calculating the variance of the phase derivative is as follows: (14) In the formula, Represents pixels The phase derivative variance This represents the partial derivative of the wrapping phase in the x-direction; This represents the partial derivative of the wrapping phase in the y-direction; express exist The average value within the window; express exist The average value within the window; where, , , This indicates that the calculation result is restricted to a specific package operator; The formula for calculating the maximum phase gradient is as follows: （15） In the formula, Represents pixels The maximum phase gradient.

6. The tooth surface phase unwrapping method based on adaptive unscented information filtering according to claim 2, characterized in that, Step 5, the method for setting the threshold, specifically includes the following sub-steps: Step 5.1: Calculate the residual points in the package phase, construct a gray value histogram, and record the location of the residual points and the corresponding gray values; Step 5.2: Count the number of residual points corresponding to each gray value, and establish the relationship between the percentage of residual points and the gray value; Step 5.3: From the percentage of residual points Starting from 0, according to The amplitude increases incrementally, and the grayscale value is recorded. ; Step 5.4: Stop increasing the threshold when the following formula condition is met. The gray value corresponding to the residual percentage at this point is the threshold: （16） In the formula, k Indicates the number of times; The non-interferometric image is divided according to the following conditions: the region composed of pixels with gray values ​​greater than or equal to the threshold is the reliable region, and the region composed of pixels with gray values ​​less than the threshold is the unreliable region; at the same time, since the non-interferometric image corresponds one-to-one with the pixels in the wrapped phase image, the wrapped phase image is also divided into reliable and unreliable regions.

7. The tooth surface phase unwrapping method based on adaptive unscented information filtering according to claim 2, characterized in that, Step 7 specifically includes: Unpack the pixels to be unpacked The initial values ​​of the state variables and state covariance matrix are estimated by incorporating them into the unscented information filter, as shown below: （17） In the formula, Represents state variables Initial value; Represents the state covariance matrix Initial value; Indicates pixels that have been unwrapped. The weights; This indicates that the unwrapped pixels have been identified. State variables; This indicates that the unwrapped pixels have been identified. The covariance matrix; Represents pixels The corresponding covariance matrix; Represents pixels Relative to pixels The local phase gradient estimate; Indicates the signal-to-noise ratio; Represents pixels Signal-to-noise ratio; Represents pixels The set of unwrapped pixels within the neighborhood; Initialize the state variables Initial values ​​of the state covariance matrix The corresponding Sigma point set is obtained by performing an unscented transformation, as follows: （18） In the formula, Represents the Sigma point. d This represents the dimension of the state variables. The scaling parameter of the Sigma point is represented by, where and For tuning factor, Take 0, ; The weights corresponding to each Sigma point are calculated as follows: （19） In the formula, They represent The weights of the state variables and the state covariance matrix; Indicates respectively The weights of the state variables and the state covariance matrix; Indicates the regulating factor. Take 2; Based on the calculated Sigma point set, for each pixel... The prediction and estimation of the state variables, as well as the corresponding prediction and estimation of the state covariance matrix, are calculated as follows: （20） In the formula, This represents the predicted estimate of the state variable. This represents the predicted estimate of the state covariance matrix; Based on the definition of information filtering, calculate the information matrix. and information vector , means as follows: （21） For pixels To predict the measurement vector, substitute the Sigma point set into the measurement equation and calculate the measurement vector covariance matrix and cross-covariance matrix. The calculation formulas are as follows: （22） In the formula, This represents the measurement value of each point calculated using the Sigma point set; This indicates that pixel points are estimated using the Sigma point set. The measured value, i.e. the predicted measured value; For the measurement equation, ; The cross-covariance matrix representing the measurement vectors; This represents the covariance matrix of the measurement vector.

8. The tooth surface phase unwrapping method based on adaptive unscented information filtering according to claim 7, characterized in that, Step 8 specifically includes: The innovation is the difference between the actual measurement and the predicted measurement, expressed as follows: （23） In the formula, Indicates new information; Indicates the measured value; Indicates the predicted measurement value; Based on the Mahalanobis distance and the statistics of the information obtained during the filtering process, and given that they follow a chi-square distribution, it is expressed as follows: （24） In the formula, Indicates the transposition of new information; Indicates the measurement noise covariance; A statistical measure representing new information; For the dimension of the new information, This indicates that the distribution follows a chi-square distribution; The threshold is determined by the chi-square distribution of the new information, and the measurement noise adjustment coefficient is calculated as follows: （25） （26） In the formula, Indicates the measurement noise adjustment factor; Represents pixels The amount of new information statistics; This represents the updated measurement noise covariance; , , Represents pixels Signal-to-noise ratio at the location.

9. The tooth surface phase unwrapping method based on adaptive unscented information filtering according to claim 8, characterized in that, Step 9 specifically includes: The information contribution is calculated using the following formula: （27） In the formula, Observational contribution items representing information state Represents the observation contribution term of the information matrix; For pixels Phase expansion is performed, as follows: （28） In the formula, , , Representing pixels The updated information state, the updated information matrix, and the final covariance matrix. Represents pixels The unfolded phase.

Images