Motion error compensation method for fringe projection phase shift profilometry

Abstract

The application discloses a motion error compensation method for fringe projection phase shift profilometry, and can realize high-precision phase recovery of an object moving in any direction. Firstly, the first fringe image of a three-step phase shift method is taken as a reference to track the motion of the measured object on the image in the second and third fringe images, and the motion compensation is performed on the second and third fringe images according to the estimated pixel translation amount, so that the object is in the same position in the three images, and the phase shift error caused by the motion becomes uniform as a whole. Then, the unknown phase shift extraction algorithm is used to estimate the new phase shift value of the compensated second and third fringe images, and the wrapped phase is calculated in combination with the three fringe images and the estimated new phase shift amount. Compared with the traditional phase shift profilometry, the application can greatly reduce the phase error caused by the motion, and realize high-precision phase recovery of the object moving in any direction.

Description

Technical Field

[0001] This invention belongs to the field of optical measurement technology, specifically a motion error compensation method for fringe projection phase-shifting contouring. Background Technology

[0002] With the advancement of optoelectronic technology, optical 3D measurement technology has also developed rapidly. Fringe projection technology, as a type of optical 3D measurement technology, features high precision, high spatial resolution, and fast measurement speed. Currently, the main fringe analysis methods include Fourier profilometry (see "Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry," by M Takeda et al.) and phase-shifting profilometry (see "Temporal phase unwrapping algorithms for fringe projection profilometry: A comparative review," by C Zuo et al.). Fourier profilometry is a single-frame 3D imaging technique that extracts fringe phase information from a single grating fringe pattern through spatial frequency filtering. Therefore, this method has high measurement efficiency and is insensitive to motion, making it suitable for 3D measurement of dynamic scenes. However, due to spectral aliasing and spectral leakage, Fourier profilometry has low fidelity in preserving contour details, resulting in low measurement accuracy and making it difficult to use for measuring objects with complex surface topography. Although windowed Fourier transform profilometry ("Two-dimensional windowed Fourier transform for fringe pattern analysis: principles, applications and implementations", by Q Kemao et al.) and wavelet transform profilometry ("Spatial carrier-fringe pattern analysis by means of wavelet transform: wavelet transform profilometry", by Zhong J et al.) have been proposed to improve phase extraction accuracy, this means higher computational costs.

[0003] Phase-shifting profilometry uses multiple consecutive projections of grating fringe patterns for phase demodulation. Generally, the phase demodulation accuracy of phase-shifting methods is much higher than that of Fourier methods, and the accuracy increases with the number of phase shifts. However, this means reduced measurement efficiency and increased sensitivity to motion. Phase-shifting profilometry requires ensuring that the position of the measured object remains unchanged during multi-frame projection. If the object moves during the measurement process, motion artifacts will occur, leading to measurement errors. Therefore, phase-shifting profilometry is only suitable for static or quasi-static measurements.

[0004] Currently, neither of the two mainstream fringe analysis methods can achieve high-precision 3D measurement of dynamic scenes. While Fourier profilometry can achieve single-frame measurement, its accuracy is low; and while phase-shifting profilometry offers high accuracy, it produces severe motion artifacts when measuring moving objects. Therefore, how to remove the phase error caused by object motion in phase-shifting profilometry and achieve motion error compensation is a problem that urgently needs to be solved. Summary of the Invention

[0005] This invention proposes a motion error compensation method for stripe projection phase-shifting contouring.

[0006] The technical solution to achieve the purpose of this invention is: a motion error compensation method for stripe projection phase-shifting contouring, the specific steps of which are as follows:

[0007] Step 1: Project a standard three-step phase-shift map with sinusoidal light intensity onto the object under test, and simultaneously capture three phase-shift fringe maps modulated by the surface of the object under test. Using the first fringe map as a reference, estimate the overall movement of the object in the second and third fringe maps to obtain the pixel translation amount.

[0008] Motion compensation is performed on the second and third stripe patterns based on the pixel translation to obtain a corrected image, so that the objects in the three stripe patterns are in the same position in the image.

[0009] Step 2: Use the random phase shift extraction algorithm to estimate the phase shift of the second and third fringe patterns relative to the first fringe pattern after motion compensation;

[0010] Step 3: Combine the corrected image and phase shift, and use the least squares method to calculate the wrap-around phase map.

[0011] Preferably, the phase shifts of the three phase-shifted fringe patterns modulated by the surface of the object under test are 0, 0, 0.

[0012] Preferably, the specific method for estimating the overall movement of the object in the second and third fringe images, using the first fringe image as a reference, to obtain the pixel translation amount is as follows:

[0013] Perform Fourier transforms on the three stripe patterns respectively to obtain the corresponding spectrum diagrams;

[0014] The normalized power spectra of the second and third images and the first image were calculated respectively. The normalized power spectra were then windowed using a rectangular window to obtain the windowed phase correlation matrix.

[0015] Principal component analysis based on singular value decomposition was used to reduce the dimensionality of the windowed phase correlation matrix.

[0016] Take the phase of the left and right principal singular vectors and unwrap them. Perform least squares fitting on the unwrapped results and obtain the pixel translation amount based on the slope.

[0017] Preferably, the first stripe pattern I1 is represented as:

[0018] I1(r)=A(r)+B(r)cos(2πf0r+φ h (r))

[0019] The nth stripe pattern I n Represented as:

[0020] I n (r)=A(r+r n )+B(r+r n cos[2πf0r+φ h (r)+2π(n-1) / N+ε n (r)]

[0021] Where n = 2, 3, N = 3, r = (x, y) T r represents the camera pixel coordinates. n =(x n ,y n ) T x represents the pixel translation caused by the motion of an object in the image. n and y n r n The components in the x and y directions, f0 = (f x ,f y f represents the fringe frequency on the camera's image plane. x and f y Let f0 be the x and y components of f0, respectively, and let A and B represent the background light intensity and fringe modulation degree, respectively. h ε represents the phase modulation caused by changes in the height of the measured object. n This represents the phase shift error caused by motion.

[0022] Preferably, Fourier transforms are performed on the three fringe patterns respectively to obtain the corresponding spectrum diagrams:

[0023]

[0024]

[0025] Where i is the imaginary unit, and f represents the frequency. Representing I1 and I respectively n Spectrum diagrams of A and B This represents the convolution operation. Indicates Fourier transform;

[0026] Calculate the normalized power spectrum of the spectrograms of the second and third images and the spectrogram of the first image, respectively:

[0027] Where, ψ n (f) represents the phase correlation matrix after windowing, * indicates taking the complex conjugate, and Rect(f) indicates applying a rectangular window to the zero frequency. It is a rank-one matrix;

[0028] Principal component analysis based on singular value decomposition was used to analyze ψ n (f) Perform data dimensionality reduction, ψ n The dimension of (f) is (H,W), and H≤W;

[0029]

[0030] Where SVD[] represents principal component analysis operation, u1, u2, ... u H It is a left singular vector, v1, v2, ... v H It is a right singular vector, σ1,σ2,…σ H It is a singular value;

[0031] Take the phase of the left and right principal singular vectors u1 and v1 respectively, unwrap them, and perform least-squares fitting on the results. The slope of the fitting results is obtained as k. y and k x , through k y and k x Receive pixel translation amount r n =(x n ,y n ) T

[0032]

[0033]

[0034] Preferably, based on the estimated pixel translation amount r n For n = 2 and 3, perform motion compensation on the second and third striped images to ensure the object is in the same position in all three images:

[0035] I'n (r)=I n (rr n ).

[0036] Preferably, the specific method for estimating the phase shift of the second and third fringe patterns relative to the first fringe pattern using the random phase shift extraction algorithm in step 2 is as follows:

[0037] After image motion compensation, the first, second, and third stripe patterns are I'1, I'2, and I'3, respectively, with corresponding phase shifts of δ1, δ2, and δ3.

[0038] Since the first fringe pattern I1 is used as a reference, no motion compensation is performed, I'1 = I1, and the phase shift satisfies 0 = δ1 < δ2 < δ3 < 2π.

[0039] The difference d between I'1 and I'2 12 for:

[0040] d 12 =I'1-I'2

[0041] Define d 12 variance for:

[0042]

[0043] Where U represents the effective area of ​​the stripe pattern, M represents the number of pixels in the effective area, and (x,y) represents the number of pixels in the stripe pattern.

[0044] The difference d between I'1 and I'3 13 variance The difference d between I'2 and I'3 23 variance

[0045]

[0046]

[0047] The phase shift δ2 of the second fringe pattern I'2 and the phase shift δ3 of the third fringe pattern I3' are determined by statistical properties:

[0048]

[0049]

[0050] Preferably, the wrap-around phase map is calculated using the least squares method by combining the corrected image and the phase shift:

[0051]

[0052]

[0053] Compared with the prior art, the significant advantages of this invention are: it can greatly reduce the phase error caused by motion and achieve high-precision phase demodulation of objects moving in any direction.

[0054] The present invention will now be described in further detail with reference to the accompanying drawings. Attached Figure Description

[0055] Figure 1 This is a flowchart illustrating the motion error compensation method for fringe projection phase-shifting profilometry. (a) shows the original three-step phase-shifting fringe image in a moving scene, and (b) shows... and (c) is the phase correlation matrix between the phases, (d) is the principal component analysis result of (c), (e) represents the least squares fitting of the principal singular vectors, (f) is I2 after image compensation, (g) represents the trigonometric relationship between the estimated phase shift and variance, (h) is the wrapped phase map after motion compensation, (i) is the 3D reconstruction result after motion compensation, and (j) is the 3D reconstruction result before motion compensation.

[0056] Figure 2 The results show the three-dimensional measurement of David in a motion scene, where (a) and (b) are the three-dimensional reconstruction results of traditional phase-shifting profilometry and the method proposed in this invention, respectively.

[0057] Figure 3 The results are three-dimensional measurements of a standard ceramic sphere in a motion scenario. (a) and (b) are the three-dimensional reconstruction results of the traditional phase-shifting profilometry and the method proposed in this invention, respectively. (c) and (d) are the errors between the sphere measured by the traditional phase-shifting profilometry and the fitted sphere by the method proposed in this invention, respectively. (e) and (f) are histograms of (c) and (d), respectively. (g) is a cross-sectional view of the straight line drawn in (a) and (b). Detailed Implementation

[0058] A motion error compensation method for fringe projection phase-shifting profilometry is proposed, which can achieve high-precision phase demodulation of objects moving in any direction. Combined with... Figure 1 The specific implementation steps of the method proposed in this invention are as follows:

[0059] Step 1: Project a standard three-step phase-shift pattern with a sinusoidal light intensity distribution onto the object under test using a projector. Simultaneously capture the pattern using an industrial camera to obtain three phase-shift fringe patterns modulated by the surface of the object under test, with phase shifts of 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 10 ... Using the first phase shift fringe pattern I1 as a reference, estimate the second and third phase shift fringe patterns I. n The overall movement of the object in the image (n=2,3) is obtained by the pixel translation amount r.n (n = 2, 3). According to r n For stripe pattern I n The corrected image I is obtained by performing translation compensation. n (n=2,3) This method places the object in the same position in the three fringe patterns, changing the non-uniform phase shift error caused by motion to a uniform one. The specific process is as follows:

[0060] The first stripe pattern I1 can be represented as:

[0061] I1(r)=A(r)+B(r)cos(2πf0r+φ h (r))

[0062] The nth stripe pattern I n (n = 2, 3) can be represented as:

[0063] I n (r)=A(r+r n )+B(r+r n cos[2πf0r+φ h (r)+2π(n-1) / N+ε n (r)]

[0064] Where N is 3, and r = (x, y) T r represents the camera pixel coordinates. n =(x n ,y n ) T x represents the pixel translation caused by the motion of an object in the image. n and y n r n The components in the x and y directions, f0 = (f x ,f y f represents the fringe frequency on the camera's image plane. x and f y Let f0 be the x and y components of f0, respectively, and let A and B represent the background light intensity and fringe modulation degree, respectively. h ε represents the phase modulation caused by changes in the height of the measured object. n This represents the phase shift error caused by motion.

[0065] For the first fringe pattern I1 and the nth fringe pattern I1 respectively n Perform a Fourier transform to obtain the corresponding spectrum:

[0066]

[0067]

[0068] Where i is the imaginary unit, and f represents the frequency. Representing I1 and I respectively n Spectrum diagrams of A and B This represents the convolution operation. This represents the Fourier transform.

[0069] To extract the pixel translation amount r n phase shift term Simultaneously, the influence of fringe modulation is eliminated, the normalized power spectrum is calculated, and a rectangular window is used to window the normalized power spectrum:

[0070]

[0071] Where, ψ n (f) represents the phase correlation matrix after windowing, * indicates taking the complex conjugate, and Rect(f) indicates applying a rectangular window to the zero frequency. Because It is a rank-one matrix, used to extract ψ n (f) Principal components, excluding the interference of high-dimensional noise, are analyzed using principal component analysis based on singular value decomposition for ψ. n (f) Perform data dimensionality reduction. Assume ψ n The dimensions of (f) are (H,W), and H≤W.

[0072]

[0073] Where SVD[] represents principal component analysis operation, u1, u2, ... u H It is a left singular vector, v1, v2, ... v H It is a right singular vector, σ1,σ2,…σ H It is a singular value.

[0074] Take the phase of the left and right principal singular vectors u1 and v1 respectively, unwrap them, perform least squares fitting on the results, and obtain k by taking the slope of the fitting results. y and k x , through k y and k x The pixel translation amount r can then be obtained. n =(x n ,y n ) T :

[0075]

[0076]

[0077] Based on the estimated pixel translation amount r n (n=2,3), perform motion compensation on the 2nd and 3rd striped images to make the object appear in the same position in all three images:

[0078] I' n (r)=I n (rr n )

[0079] Step 2: Use the random phase shift extraction algorithm to estimate the compensated second and third stripe images I' n (n=2,3) Phase shift δ relative to the first fringe pattern I1 n The specific process is as follows:

[0080] Assume that the 1st, 2nd, and 3rd fringe patterns after image motion compensation are I'1, I'2, and I'3, respectively, with corresponding phase shifts of δ1, δ2, and δ3. Since the 1st fringe pattern I1 is used as a reference and no motion compensation is performed, I'1 = I1, and the phase shifts satisfy 0 = δ1 < δ2 < δ3 < 2π.

[0081] The difference d between I'1 and I'2 12 for:

[0082] d 12 =I'1-I'2

[0083] Define d 12 variance for:

[0084]

[0085] Where U represents the effective region of the stripe pattern, M represents the number of pixels in the effective region, and (x,y) are the number of pixels in the stripe pattern. Similarly, the difference d between I'1 and I'3 can be obtained. 13 variance The difference d between I'2 and I3' 23 variance

[0086]

[0087] The phase shift δ2 of the second fringe pattern I'2 and the phase shift δ3 of the third fringe pattern I'3 can be obtained from the statistical properties:

[0088]

[0089]

[0090] Step 3, combine the three striped images I' from Step 1 after image motion compensation. i (i = 1, 2, 3) and the estimated phase shift δ in step 2 i (i = 1, 2, 3), the wrapping phase map φ is calculated using the least squares method. The specific process is as follows:

[0091] Using the least squares method, we can obtain:

[0092]

[0093]

[0094] Example

[0095] To verify the effectiveness of this method, a fringe projection 3D measurement system was built. This system includes one DLP projector (LightCrafter 4500Pro) and three monochrome cameras (Basler acA640-750um). The projector has a resolution of 912×1140, and the monochrome cameras have a resolution of 640×480. Each camera is equipped with a 12mm lens. The working distance between the system and the object being measured is approximately 1m. The projection and imaging speed is set to 100 frames per second. The projector projects a standard three-step phase-shift map, and the number of periods of the projected fringe pattern is set to 64. The system uses a stereo phase unwrapping algorithm based on multi-view geometric constraints (reference “High-precision real-time 3D shape measurement based on a quad-camera system,” author Tao T) for absolute phase unwrapping.

[0096] The effectiveness of the method described in this invention was first verified by measuring the moving plaster cast of David. The plaster cast of David was moved at a speed of 80 mm / s using a displacement stage. The three-dimensional reconstruction results obtained using conventional phase-shifting profilometry and the method described in this invention are as follows: Figure 2 As shown in (a) and (b), it can be seen that the 3D surface reconstructed using the traditional phase-shifting profilometry has obvious motion artifacts, while the 3D surface obtained using the method described in this invention removes motion artifacts and has a higher quality reconstruction result.

[0097] The second motion scenario involves a standard ceramic ball moving at approximately 80 mm / s. The measurement results are as follows: Figure 3As shown in the figure, (a) and (b) are the 3D reconstruction results obtained by the traditional phase-shifting profilometry and the method described in this invention, respectively; (c) and (d) represent the errors between the sphere measured by the traditional phase-shifting profilometry and the fitted sphere by the method described in this invention, respectively; (e) and (f) are the histograms of (c) and (d), respectively; and (g) is a cross-sectional view of the straight line drawn in (a) and (b). The results show that the sphere measured by the traditional phase-shifting profilometry has obvious ripples and a large sphere fitting error, with an accuracy of 0.33268 mm; the sphere measured using the method described in this invention has no motion ripples and a very small sphere fitting error, with an accuracy of 0.07688 mm. The measurement results of the above scene demonstrate that the method proposed in this invention can suppress motion artifacts and achieve high-precision phase demodulation and 3D measurement of moving objects.

Claims

1. A method for motion error compensation in striped projection phase-shifting contouring, characterized in that, The specific steps are as follows: Step 1: Project a standard three-step phase-shift map with a sinusoidal light intensity distribution onto the object under test. Simultaneously capture three phase-shift fringe maps modulated by the surface of the object under test. Using the first fringe map as a reference, estimate the overall movement of the object in the second and third fringe maps to obtain the pixel translation. The specific method is as follows: Perform Fourier transforms on the three stripe patterns respectively to obtain the corresponding spectrum diagrams; The normalized power spectra of the second and third images and the first image are calculated respectively. The normalized power spectra are then windowed using a rectangular window to obtain the windowed phase correlation matrix. Principal component analysis based on singular value decomposition was used to reduce the dimensionality of the windowed phase correlation matrix. Take the phase of the left and right principal singular vectors and unwrap them. Perform least squares fitting on the unwrapped results and obtain the pixel translation amount based on the slope. Motion compensation is performed on the second and third stripe patterns based on the pixel translation to obtain a corrected image, so that the objects in the three stripe patterns are in the same position in the image. Step 2: Use the random phase shift extraction algorithm to estimate the phase shift of the second and third fringe patterns relative to the first fringe pattern after motion compensation; Step 3: Combine the corrected image and phase shift, and use the least squares method to calculate the wrap-around phase map.

2. The motion error compensation method for stripe projection phase-shifting contouring as described in claim 1, characterized in that, The phase shifts of the three phase-shifted fringe patterns modulated by the surface of the object under test are respectively... , , .

3. The motion error compensation method for stripe projection phase-shifting contouring as described in claim 1, characterized in that, The first striped pattern Represented as: No. striped pattern Represented as: in, , It is 3. Represents the camera pixel coordinates, This represents the amount of pixel translation caused by the motion of an object in the image. and They are respectively exist and Components in direction, This represents the frequency of the fringes on the camera's image plane. and They are respectively exist and Components in direction, , These represent the background light intensity and the stripe tone, respectively. This indicates the phase modulation caused by changes in the height of the measured object. This represents the phase shift error caused by motion.

4. The motion error compensation method for stripe projection phase-shifting contouring as described in claim 3, characterized in that, Perform Fourier transforms on the three fringe patterns respectively to obtain the corresponding spectrum diagrams: in, It is the imaginary unit. Represents frequency, , , , They represent , , , The spectrum diagram, This represents the convolution operation. Indicates Fourier transform; Calculate the normalized power spectrum of the spectrograms of the second and third images and the spectrogram of the first image, respectively: in, This represents the phase correlation matrix after windowing. Indicates taking the conjugate. This indicates adding a rectangular window to the zero frequency. It is a rank-one matrix; Principal component analysis based on singular value decomposition was used to... Perform data dimensionality reduction. The size is ,and ; in, This represents principal component analysis. It is a left singular vector. It is a right singular vector. It is a singular value; For the left and right principal singular vectors respectively and Take the phase and unwrap it, then perform a least-squares fit on the result. The slopes of the fitted results are obtained respectively. and ,pass and Receive pixel translation amount 。 5. The motion error compensation method for stripe projection phase-shifting contouring as described in claim 1, characterized in that, Based on the estimated pixel translation , Motion compensation is applied to the second and third striped images to ensure the object is in the same position across all three images. 。 6. The motion error compensation method for stripe projection phase-shifting contouring as described in claim 1, characterized in that, Step 2, which uses the random phase shift extraction algorithm to estimate the phase shift of the second and third fringe patterns relative to the first fringe pattern after motion compensation, is as follows: The first, second, and third stripe patterns after image motion compensation are respectively , , The corresponding phase shifts are respectively , , ; Due to the first stripe pattern For reference only, no exercise compensation is performed. And the phase shift satisfies and difference for: definition variance for: in, This indicates the effective area of ​​the stripe pattern. Indicates the number of pixels in the valid area. The pixels of the stripe pattern; and difference variance and and difference variance : The second stripe pattern was determined by statistical properties. phase shift And the third striped pattern phase shift : 。 7. The motion error compensation method for stripe projection phase-shifting contouring as described in claim 6, characterized in that, By combining the corrected image and the phase shift, the wrap-around phase map is calculated using the least squares method: 。

Images