Phase deflection measurement gradient field three-dimensional surface shape reconstruction method

Abstract

The invention discloses a phase deflection measurement gradient field three-dimensional surface shape reconstruction method. Aiming at the problems of low precision and poor stability when a traditional gradient integral or Poisson reconstruction algorithm is used for processing a missing measurement area, noise data and large-scale data, the invention provides a unified adaptive solution framework. The core idea of the method is that on the basis of inputting gradient field data and mask moments, a linear equation set containing gradient constraints and thin plate regularization constraints at the same time is established, and global optimal reconstruction is achieved in a least square mode. When the number of effective data points is small, direct solution is adopted; when the data scale is large, a rapid iteration solving strategy combining a precondition conjugate gradient method (PCG) and an incomplete Cholesky (IC) precondition device is adaptively switched, smooth compensation of a missing measurement area can be guaranteed, the detail precision of the curved surface can be maintained, the calculation efficiency is greatly improved, and the method is suitable for large-scale measurement. The method is high in precision, rapid in convergence and good in robustness in measurement of small workpieces and large free-form surfaces.

Description

Technical Field

[0001] This invention mainly relates to the field of three-dimensional topography reconstruction and surface measurement technology, specifically to an adaptive three-dimensional surface shape recovery method and system based on gradient field and thin plate regularization. Background Technology

[0002] The gradient field and thin-plate regularization adaptive 3D surface reconstruction method, as a high-precision surface reconstruction technique, is widely used in optical measurement, precision manufacturing, reverse modeling, 3D inspection, and computer vision. This type of method establishes corresponding gradient constraint equations by measuring the gradient field or normal vector field of the object surface, and combines this with a thin-plate regularized smoothing model to globally optimize and reconstruct the surface topography. By introducing a regularization term during the optimization process, the influence of noise can be effectively suppressed and the continuous topography of the missing measurement areas can be restored, achieving high-precision surface fitting and smooth reconstruction. However, traditional methods often face problems such as high computational complexity, large memory consumption, and insufficient reconstruction accuracy in the missing measurement areas when processing large-scale data or high-noise scenes, resulting in low algorithm convergence efficiency and distortion or discontinuity in the reconstruction results.

[0003] Traditional methods for 3D surface reconstruction mainly include direct integration, the Poisson equation method, least squares fitting, and their improved algorithms. Direct integration obtains the height value by integrating the gradient field point by point. For example, Frankot and Chellappa et al. proposed projecting the non-integrable gradient field into an integrable subspace in the Fourier domain, and then recovering the height field using frequency domain integration. However, this method is sensitive to noise and boundary conditions, and projection based on global basis functions is prone to artifacts or integration drift when there are large areas of missing measurements or non-periodic boundaries.

[0004] The Poisson equation method transforms the gradient field reconstruction problem into solving a two-dimensional Poisson equation. Represented by Pérez et al.'s Poisson image editing and Simchony et al.'s reduction of gradient domain reconstruction to solving a two-dimensional Poisson equation, the height is recovered from the gradient field by solving (or least squares) the equation. However, in cases of missing data, uncertain boundary conditions, or incompletely integrable gradient fields, integral drift can occur. Direct solutions are computationally and memory-intensive for large-scale meshes, necessitating optimization using sparse / iterative methods.

[0005] Least squares fitting solves for the globally optimal surface by constructing an overdetermined system of equations. Harker and O'Leary et al. systematically proposed a regularized reconstruction framework based on matrix algebra and put forward an efficient algorithm that can be scaled to large pixel scales. Although the theory is sound, preconditioners and numerical stabilization strategies are still needed for extremely large-scale or column-dependent problems; otherwise, the solution will still be time-consuming or unstable. Summary of the Invention

[0006] The purpose of this invention is to provide a high-precision and robust three-dimensional surface shape recovery method. By combining gradient measurement constraints and thin plate regularization constraints, a unified least squares optimization problem is constructed, which effectively suppresses noise interference, handles data missing regions, and adaptively selects the optimal solution strategy according to the data scale, thereby achieving fast, stable and accurate reconstruction of the three-dimensional surface shape of an object.

[0007] The objective of this invention is achieved as follows: First, according to the method described in claim 1, data preprocessing and index allocation are performed.

[0008] Step 1-1: Input gradient measurement matrix G x G y With the mask matrix mask, at this time G x (i,j), G y (i,j) represent the gradient measurements of pixel (i,j) in the x and y directions, respectively. mask(i,j)=1 indicates that the point is a valid measurement point, and 0 indicates that it is invalid or missing measurement point.

[0009] Steps 1-2: After dividing the valid / invalid regions, assign a global index (k=1,2,…,N) to all points (i,j) that satisfy mask(i,j)=1, and establish a mapping table from k to (i,j) to map the solution z in vector form back to two-dimensional pixel coordinates. At the same time, count the total number of valid points N.

[0010] Secondly, according to the method described in claim 2, the gradient constraint matrix A is constructed. g (Differential constraints in the x and y directions).

[0011] Step 2-1: For each valid point (i,j), a difference equation can be established only if the adjacent pixels are also valid points: ;

[0012] Where z(i,j) represents the height of pixel (i,j), this data is to be determined, G x (i,j), G y (i,j) corresponds to the matrix input in step 1-1 of claim 1, and represents the gradient measurement values ​​of the point in the x and y directions, respectively.

[0013] Step 2-2: Rearrange all difference equations into a sparse linear system: ; Where, matrix A gThis is a sparse gradient constraint matrix, where each row corresponds to a difference constraint equation. Each row contains only two non-zero terms, for example, +1 at column z(i+1,j) and -1 at column z(i,j). z is the variable to be solved, with a length equal to the height vector N described in steps 1-2 of claim 1. g G is a column vector of gradient observations. x or G y The corresponding entry.

[0014] Then, according to claim 3, the thin plate regularization smoothing matrix A is constructed. r .

[0015] Step 3-1: For a non-boundary valid point (i,j) (where all four of its neighbors are valid), establish a five-point second-order difference constraint: ; Among them, the four neighborhood terms z(i±1,j) and z(i,j±1) represent the pixel heights adjacent to the center point, the coefficient -4 represents the second-order bending penalty of the center point, and the right side of the equation is 0, indicating that the expected second-order difference is close to 0, thus suppressing high-frequency oscillations.

[0016] Step 3-2: Gather all smoothing constraints and construct a sparse matrix: ; Among them, A r This represents a sparse smooth constraint matrix, where each row is a five-point template with a center coefficient of -4 and four neighbor coefficients of +1. r It is usually a zero vector, or the corresponding value when it contains boundary / prior values.

[0017] Next, according to claims 1 and 4, the constraints are combined and a unified optimization problem is constructed: Step 4-1: Set the regularization parameter to λ>0, and construct the extended matrix and vector: , ; in, This serves to balance the weights of the two items at the matrix level. In this case, the 0 in the lower half corresponds to b in step 3-2 of claim 3. r If b r If ≠0, replace it with That's all.

[0018] Step 4-2: To ensure that the final reconstructed surface both conforms to the measurement data and maintains overall smoothness, we combine these two requirements by establishing a least-squares problem: ; Where || represents the Euclidean norm, which means squaring and summing these errors; Az-b in the formula represents the error between the calculated value and the measured value, and min at the beginning means minimizing the error.

[0019] Step 4-3: The z in the previous step cannot be directly calculated; it needs to be transformed into a form that can be directly solved using matrix operations. Mathematically, minimizing the error is equivalent to making its derivative zero, that is, when the error is minimized, the partial derivative with respect to z is 0. Expanding and simplifying the expression, we get: ; Let A T A=M, representing the combined weight of all constraints, let A T b=r, representing the target trend derived from the observed data, so the formula can be written as: ; Solving for z yields the height distribution that simultaneously satisfies both gradient and smoothness constraints.

[0020] Subsequently, an adaptive solution method is selected and the solution is obtained. There are two methods: one is the direct solution method, and the other is the solution method that combines the preconditioned conjugate gradient method (PCG) with the incomplete Cholesky (IC) preconditioner.

[0021] In step 1-2 of claim 1, if the number of valid points N is less than N th At this time N th If the preset threshold is met, the direct solution process will be initiated; otherwise, the iterative solution process will be initiated.

[0022] First, there's the direct solution method suitable for small-scale applications, which only requires direct least squares solutions: ; (A T A) -1 This method represents the inversion of a matrix when it is numerically invertible, or it can be achieved through LU or QR decomposition. It is suitable for cases where N is small and the matrix conditions are good. Although this method involves slightly more computation, it can be completed in a short time and the results are accurate.

[0023] Secondly, there is an adaptive iterative solution method suitable for large-scale applications, namely the preconditioned conjugate gradient method (PCG), which is combined with an incomplete Cholesky (IC) preconditioner to accelerate convergence.

[0024] Step 5-1: Forming the system matrix , ; Step 5-2: Try using Cholesky decomposition or detecting the minimum eigenvalue to determine if it is symmetric positive definite (SPD). If it is SPD, perform incomplete Cholesky (IC) decomposition: ; Where L is a lower triangular sparse matrix, IC preserves sparsity to save memory, and P = LL is set. T As a preconditioner, it is an approximate matrix obtained through incomplete Cholesky decomposition (IC).

[0025] Step 5-3: We first use an easily manageable matrix P to correct the shape of the original equation, making it easier to solve using the conjugate gradient method: ; At this time, matrix P -1 M becomes easier to compute. Next, we proceed with iterative solving; the algorithm's loop can be simply summarized as follows: Step 5-3-1: Perform initialization, using z0 as the initial solution, calculate the residual s0=r-Mz0, and set a direction vector p0=P -1 s0.

[0026] Step 5-3-2: Perform iterations, first calculate the step size α. k Then update the solution vector: ; Then recalculate the current error: ; Finally, determine whether convergence has occurred, i.e., whether the error is small enough: ; Where ε is the convergence threshold. If the condition is met, the result is sufficiently accurate, and the calculation stops; if convergence fails, the direction vector p is updated. k+1 And continue iterating.

[0027] Sometimes matrix M is not a standard symmetric positive definite (SPD) matrix, which can cause the conjugate gradient method to malfunction. In such cases, this invention automatically corrects the matrix by adding a small diagonal term: ; This can help the system regain stability.

[0028] Finally, according to claim 1, the output is mapped and post-processed.

[0029] Step 6-1: After the previous optimization and solution, a vector-like result was obtained: ; It contains the height values ​​of all valid measurement points, that is, the height information of each valid pixel on the three-dimensional surface. At this point, z needs to be placed back to the corresponding two-dimensional coordinate position in order to recover the complete three-dimensional shape.

[0030] In the data preprocessing stage of steps 1-2 of claim 1, a unique global index k has been assigned to each valid point, and a mapping table between index k and coordinates (i,j) has been established. This mapping table is then used to map the k-th element z in the solution vector. k Assigning the pixel coordinates back to their corresponding (i,j) results in the two-dimensional height distribution matrix: ; Step 6-2: In the measurement data, there are some invalid or missing areas, such as surface reflection, occlusion, or detection blind spots. The pixels in these areas are empty when mask=0, and these holes need to be smoothed and filled.

[0031] This invention employs two smooth filling strategies: the first is a neighborhood interpolation method, which estimates the missing measurement points based on the average height of adjacent valid points; the second is a regularization-based smooth extension, which utilizes the thin plate regularization matrix A. r Apply a smoothing constraint to the missing measurement area so that the reconstructed surface transitions naturally at the edges.

[0032] After obtaining the complete height distribution, the results need to undergo certain geometric and numerical post-processing to ensure their accuracy and comparability. Possible methods include plane fitting and deslant correction, scale normalization, or comparison with contact measurements and calculation of RMS error.

[0033] Step 6-3: Finally, the z-series is processed. map (i,j) can be exported as point cloud data (XYZ format), mesh surface model (STL, PLY, etc.) or 2D grayscale height map, and generate visualizations, error reports or data formats for subsequent processing and inspection. Attached Figure Description

[0034] Figure 1 This is the overall algorithm flowchart of the present invention.

[0035] Figure 2 This is a framework diagram for reconstructing the three-dimensional surface shape of a test object by measuring the gradient field of a phase deflection measurement.

[0036] Figure 3 It is a three-dimensional schematic diagram of the real curved surface of the test piece.

[0037] Figure 4 It is a theoretical phase distribution map obtained by linear mapping of a three-dimensional surface.

[0038] Figure 5 It is a simulated fringe pattern generated from the theoretical phase.

[0039] Figure 6 It is a wrap-around phase distribution map obtained by simulating the phase of the stripes using inverse cosine.

[0040] Figure 7 It is a continuous phase distribution map obtained after phase recovery using a two-dimensional phase unpacking algorithm.

[0041] Figure 8 It is the effective measurement region mask (binary mask image) read from the image.

[0042] Figure 9 This is a gradient field distribution map with holes calculated under mask constraints.

[0043] Figure 10 This is a three-dimensional surface morphology result image reconstructed based on the method of this invention.

[0044] Figure 11 This is a residual error distribution diagram after planar correction. Detailed Implementation This embodiment takes the surface morphology of the test piece as the research object and uses an imaging model based on the phase deflection principle to simulate the actual measurement process. This includes five steps: fringe pattern generation, phase extraction, phase unpacking, gradient calculation, and three-dimensional surface reconstruction. The reconstruction algorithm based on gradient constraints and thin-plate regularization proposed in this invention is used to restore the surface, verifying the stability, accuracy, and robustness of the method. The specific implementation process is as follows: Figure 2 As shown.

[0045] S1: Mesh Definition and Realistic Surface Construction like Figure 3 As shown, a regular mesh of size is first established, with a mesh spacing of dx = dy = 1. The true surface shape is defined as: ; This surface serves as a simulation of the actual height distribution on the workpiece surface, and is used for the subsequent generation of fringe patterns and gradients.

[0046] S2: Simulated phase generation Construct the deflection phase based on the actual surface shape: ; Where k x k y Let k be the system sensitivity constant. In this embodiment, k is taken as... x =1,k y =1.

[0047] The generated phase is wrapped using the interval [-π, π) to obtain: ; Here, "mod" refers to the modulo operation, also called the remainder operation, which restricts the original phase ф(x,y) to the interval [-π,π). Figure 4 .

[0048] S3: Simulated stripe image generation Based on the standard fringe imaging model: ; Where, the background light intensity I0=1, and the modulation degree I m =0.8 generates a single fringe pattern to simulate the fringe acquisition process of a real optical system, such as Figure 5 .

[0049] S4: Phase Extraction and Unpacking Phase extraction, i.e., arctangent demodulation, is performed on the simulated fringe pattern to obtain the following: Figure 6 Wrap phase ф w (x,y).

[0050] Unpacked continuous phase ф u (x,y) is obtained by minimizing the gradient consistency error, and its optimization model can be expressed as: ; Among them, ф w (x,y) represents the wrapping phase, ▽ф represents the gradient of the spatially continuous phase, and wrap() is the wrapping operator that constrains the error to the interval [-π,π). For example... Figure 7 As shown, this step is equivalent to a numerical simulation of "phase unpacking" in real measurements.

[0051] S5: Gradient Calculation and Noise Simulation Based on the unpacking phase, the gradients in the x and y directions are obtained by inverting the system sensitivity: , ; To simulate real measurement noise, Gaussian white noise with a standard deviation of σ = 0.01 is superimposed on the gradient, resulting in... , In addition, such as Figure 8 As shown, a hexagonal hollow annular mask is used to divide the effective measurement area and the missing measurement area. The area with a value of 1 in the mask represents the effective measurement point, and the area with a value of 0 in the mask represents the missing measurement area.

[0052] S6: 3D Surface Reconstruction Based on Gradient and Thin Plate Regularization Within the effective region of the mask, calculate the gradient components G of the surface in the x and y directions, respectively. x G y , to obtain Figure 9The diagram shows the distribution of the gradient field with holes. We construct the gradient constraint matrix and the thin-plate regularization matrix to form a unified optimization problem: ; The regularization parameter is set to λ = 0.001. When the number of valid points exceeds a preset threshold, the preconditional conjugate gradient (PCG) algorithm is used to solve the problem. ; Where M=A T A, r=A T b, at this time P≈LL T The preconditioning matrix obtained from the incomplete Cholesky decomposition is used to improve the iterative convergence, ultimately yielding the reconstructed surface shape z. rec .

[0053] S7: Reconstruction Results and Error Analysis The reconstructed three-dimensional surface morphology is as follows Figure 10 As shown, after comparing the reconstructed results with the actual surface shape and performing planar correction, see... Figure 11 The residuals exhibit a random distribution, with the root mean square error (RMS) remaining at the order of magnitude. Even in the presence of noise and missing measurement areas, the reconstruction results remain smooth and continuous, without discontinuities or edge distortions, demonstrating the stability and robustness of the method of this invention.

Claims

1. A method for reconstructing the three-dimensional surface shape of a gradient field measured by phase deflection, comprising the following steps: Step 1: Divide the gradient measurement data into valid and invalid regions, assign a global index to the points in the valid region, supplement the two sets of fringe data obtained from the phase deflection device, calculate the phase, perform phase unrolling, and obtain gradient data in two directions; Step 2: Based on the measured gradient data, establish the corresponding gradient constraint matrix; Step 3: Construct a smoothing constraint matrix based on the thin plate regularization principle; Step 4: Establish a unified least squares optimization model based on the gradient constraint matrix and the smoothing constraint matrix; Step 5: Adaptively select the solution method based on the number of valid data points; Step 6: Map the obtained result vector back to the spatial position corresponding to the effective area, eliminate the adjustment error term, and obtain complete and continuous three-dimensional surface data to achieve smooth reconstruction of the missing area.

2. The method according to claim 1, characterized in that, The construction of the gradient constraint matrix includes differential constraints on the x and y directions; Step 2-1: For each valid pixel (i,j), where i and j represent the two-dimensional pixel indexes, assign a global index number k∈(1,2,…,N) according to the rules, where N is the total number of valid pixels, and establish an index mapping relationship: ， This index number is used to uniquely identify the position of each valid pixel in the height vector to be calculated when constructing the matrix; at this time, the vector z=[z1,z2,…,z…] N ] T , which corresponds to the set of height values ​​of all valid points, where z k =z(i,j); Step 2-2: Based on the measured gradient data G x (i,j), G y For (i,j), establish the following linear equations respectively: ， ， Here, z(i,j) represents the height value of the i-th and j-th pixels. Through the aforementioned global index mapping relationship, the height value z(i,j) of each two-dimensional position (i,j) has a unique index position k in the vector z. Therefore, the above equation can be converted into a linear algebraic form. Steps 2-3: Within the effective region, for each effective pixel (i,j), measure G based on its gradient in the direction of . x (i,j), G y (i,j) can be uniformly written as the following difference relation: ， Here, "neighbor" represents a valid pixel adjacent to (i,j) in the x or y direction, and d(i,j) is the measurement gradient in the corresponding direction multiplied by the pixel spacing, such that d(i,j) = ΔxG. x (i,j) or d(i,j)=ΔyG y (i,j); Vectorizing the above difference relationship yields a unified matrix expression: ； Where, matrix A g The gradient constraint matrix has two non-zero elements in each row, which are A and J, respectively, located in the global index column corresponding to (i,j) and its neighbors. g (r,k neighbor )=+1 and A g (r,k i,j ) = -1, while b g Let b be the column vector of gradient observations. g (r)=d(i,j), z(i,j) is the height vector to be solved.

3. The method according to claim 1, characterized in that, The thin plate regularization constraint adopts the form of a discrete Laplacian operator to constrain surface curvature and enhance smoothness; Step 3-1: For non-boundary valid points (i,j), establish second-order difference constraints: ； Among them, the four neighborhood terms z(i±1,j) and z(i,j±1) represent the pixel heights adjacent to the center point, and the coefficient -4 represents the second-order bending penalty of the center point; Step 3-2: Discretize the smoothing constraint to form a sparse matrix A r At this point, the vectorized result is a regularized matrix relation A. r z=0, and its structure satisfies the five-point template form; Step 3-3: The smoothing constraint term is added to the optimization objective function through thin plate energy regularization: ， Where λ is the regularization parameter, used to control the balance between fitting accuracy and smoothness, the first term ||A g zb g || 2 The second term, ||A|, represents the gradient fitting error. r z|| 2 This indicates a penalty for a smooth surface.

4. The method according to claim 1, characterized in that, The adaptive selection solution method includes: Step 4-1: Calculate the number of valid data points N; Step 4-2: When the number of valid points N is less than the preset threshold N th In this case, the direct solution method is used: ； in, For a combined matrix, For a combined vector, (A) T A) -1 A T This is the least squares solution matrix; Step 4-3: When the number of valid points N is greater than the preset threshold N th In this case, the iterative solution method (PCG) is used, and the solution is accelerated by the preconditional conjugate gradient method: ； Where M=A T A is the system matrix, r=A T b is the equivalent right-hand side term, and P is the preconditioning matrix used to improve convergence; in this case, P = LL. T , is the lower triangular sparse matrix obtained by performing an incomplete Cholesky decomposition on M.

5. The method according to claim 4, characterized in that, The preconditioning conjugate gradient method uses incomplete Cholesky (IC) decomposition to construct the preconditioner; Step 5-1: Determine if matrix M = A T Is A a symmetric positive definite (SPD)? Step 5-2: If it is an SPD, then perform an incomplete Cholesky decomposition: ； Where L is a lower triangular sparse matrix, L T It is its transpose matrix; Step 5-3: Set P=LL T Substitute it into the PCG algorithm as a preconditioner to accelerate iterative solution; Step 5-4: If M is not SPD, then correct it: ； Here, I is the identity matrix, used to improve the positive definiteness of the matrix, the correction term λI controls the numerical stability, and finally the decomposition and solution are re-executed to ensure that the algorithm converges stably.

Images