CN103267496A - Improved window Fourier three-dimensional measuring method based on wavelet transform - Google Patents

Abstract

The invention discloses an improved window Fourier three-dimensional measuring method based on wavelet transform. The method mainly aims at solving precise phase distribution of stripe images and obtaining three-dimensional shape information of an object through the phase distribution. The method includes the steps that black-and-white sine modulation stripe images are projected on the object to be measured and collected by a CCD, the wavelet transform is carried out on the collected and deformed stripe images line by line, a wavelet transform ridge is extracted and calculated, errors caused by a phase second derivative during ridge extraction are divided, and finally a precise window size matrix is obtained. The precise window matrix is substituted into window Fourier transform, and corresponding relative phase information of the deformed stripe images is obtained after steps such as filtering. A stripe image quality diagram is established, phases are then expanded by means of the deluge algorithm, and absolute phase distribution of the stripe images is obtained. According to the phase height conversion formula, the three-dimensional information of the object to be measured can be obtained according to the absolute phase distribution.

Description

A kind of improvement window Fourier three-dimensional measurement method based on wavelet transformation

Technical field

The invention belongs to the technical field of three-dimensional information reconstruct, based on the modulated grating projection, combined with wavelet transformed and window fourier transform method are with the three-dimensional measurement of fringe projection technology of profiling for object.

Background technology

Based on the three-dimensional measurement of optical projection, be widely used in product detection and machining control, medical field, historical relic's protection field, aerospace field, culture field etc.In numerous three-dimensional measurement technology, characteristics become prevailing three-dimensional measurement technology to the optical three-dimensional measurement technology so that its non-cpntact measurement, real-time be better etc.

The optical three-dimensional measurement technology is based on contemporary optics, melts optoelectronics, Computer Image Processing, and graphics, science such as signal processing are the modern surveying technology of one.The entity that it detects optical imagery and store as information obtains the three-dimensional measurement information needed from image.With respect to other three-dimensional information measuring method, high-speed, high-precision characteristics that optical measuring technique has since the sixties in last century, have obtained studying widely and using.Three-dimensional measurement technology based on the encode grating projection is the most representative, is most widely used.

Three-dimensional measurement based on optical grating projection projects to raster pattern the measured object surface exactly, obtains the raster image of distortion by video camera, and is determined the elevation information on testee relative reference plane by phase place and the relation of height.When optical grating projection was to the body surface, periodically the phase place of grating just was subjected to the modulation of body surface height profile, formed deformed grating, had the three-dimensional information of object in the deformed grating.The PHASE DISTRIBUTION of accurately obtaining the deformed grating image plays key effect in whole three-dimensional measurement process.

Usually three mensurations based on optical grating projection are divided into multiframe processing phase-shift method and single frames processing conversion mensuration.Fourier transform profilometry (FTP) is a kind of active single frames mensuration, and the earliest by propositions such as Takeda, numerous scholars have done deep research to this.The new application of fourier transform method in the three-dimensional measurement field that Open from This Side.But because Fourier transform (FT) is global change, do not have the ability of partial analysis, so local mistake can't be analyzed solution separately.Based on this, the window Fourier transform (WFT) with partial analysis advantage is introduced into the three-dimensional measurement field.The local signal analysis makes signal processing results more accurate, and wherein key is exactly the window selection in early stage.Window is more big, and frequency resolution is more good and spatial resolution is more poor; Window is more little, and frequency resolution is more poor and spatial resolution is more good.Therefore the selection of window has significant impact for the precision of three-dimensional measurement.According to the Heisenberg uncertainty principle, we can not provide the window size of an entirely accurate, but can constantly optimize.Therefore how choosing self-adapting window becomes a difficult point, so small echo ridge method etc. are introduced into three-dimensional measurement, has provided detailed analysis for the selection of window size, and precision and the speed of window Fourier mensuration is had raising.But since window Fourier to handle high-resolution image slower, therefore have under the situation of requirement in real-time, be necessary to reduce the computing time of window size as far as possible.

Because the phase value that fourier transform method obtains obtains absolute phase values between 0-2 π so need carry out phase unwrapping.But because the window Fourier transform is asked for its window size in the general ridge method of object height saltus step and precipitous zone and is had very mistake, error can appear in the process that causes phase place to be found the solution, this error can influence the next phase unwrapping precision of pending point when phase unwrapping, so produce continuous mistake, i.e. backguy phenomenon in the time of can causing phase unwrapping.The scanning line method speed of directly launching phase place is very fast, but robustness a little less than, very easily produce the backguy phenomenon.Dual-frequency method contains the sinusoidal grating of two kinds of frequency components by projection, through Fourier transform filtering, obtain two kinds of phase values under the frequency, with the phase value that obtains under the phase value correction high frequency that obtains under the low frequency, thus the improper value at the discontinuous place of phase place that obtains under the correction high frequency.But because dual-frequency method also is to find the solution initial phase with Fourier transform, so also exist two kinds of spectral aliasing problem between the frequency component.

Summary of the invention

Goal of the invention: at wavelet transformation when asking for window Fourier window size speed slowly and height saltus step and precipitous zone ask for the big problem of window size error, the objective of the invention is to not influence window size asks under the prerequisite of accuracy, solve hop region and ask for the inaccurate problem of phase place, a kind of improvement window Fourier three-dimensional measurement method based on wavelet transformation is provided.

Technical scheme: a kind of improvement window Fourier three-dimensional measurement method based on wavelet transformation comprises the steps:

Step 1) to the testee surface, uses the CCD camera that the testee surface is taken the black and white strip image projection, obtain a panel height degree and be c, width and be r deforming stripe image g (x, y):

G (x, y)=A (x, y)+B (x, y) cos [2 π f ₀X+ φ (x, y)] (1) formula

Wherein, (x y) is the background light distribution to A, and (x y) is the body surface reflectivity, f to B ₀Be the sine streak frequency, (x y) is relative phase distribution plan to be asked to φ, (x, y) two-dimensional coordinate of expression deforming stripe image;

Step 2), to described deforming stripe figure g (x y) does wavelet transformation line by line, obtains the approximate scale factor of deforming stripe image every bit, and detailed process is as follows:

Step 2.1), y is considered as constant y ₁, adopt one-dimensional wavelet transform to described deforming stripe image g (x, y y) ₁Row is handled, and processing procedure is: $W_{y_1}(a_1,b)={\∫}_{-\∞}^{+\∞}g(x,y_1){M^{*}}_{a,b}\left(x\right)\mathrm{dx}$ (2) formula

Wherein, y ₁The row number of representing certain delegation; A is scale factor, and span is 10 to 90, gets a value every 0.2; B is shift factor, and span is 1 to the width r of striped, gets a value every 1, and unit is pixel; Obtain through one-dimensional wavelet transform

Be the two-dimentional complex matrix of one 400 capable r row, claim y ₁The row wavelet transform matrix, a ₁It is matrix

The line label of middle element; Wherein,

M ^* _{A, b}(x) be M _{A, b}(x) conjugate function; Wavelet function M (x) according to formula (3) obtains M ^* _{A, b}(x) expression formula is suc as formula (4):

$M\left(x\right)=\frac{1}{{\left(f_b^2\mathrm{\π}\right)}^{1/4}}\exp \left(2\mathrm{\πi}{f}_{c}x\right)\exp \left(\frac{{-x}^2}{{2f}_b^2}\right)$ (3) formula

${M^{*}}_{a,b}\left(x\right)=\frac{1}{\sqrt{a}{\left(f_b^2\mathrm{\π}\right)}^{1/4}}\exp (-2\mathrm{\πi}{f}_c\frac{x-b}{a})\exp [-\frac{{(x-b)}^2}{2f_b^2{a}^2}];$ (4) formula

Wherein, f _bBe the bandwidth of wavelet function, f _cBe the centre frequency of wavelet function, i is complex unit, and x is independent variable;

Step 2.2), ask for stripe pattern g (x, approximate scale factor distribution plan f y) _A1(x, y ₁):

Calculate two-dimentional complex matrix

The modular matrix of correspondence

The search modular matrix The element of x column mean maximum, and with modular matrix

The line label assignment of x column mean greatest member give a _Max, the approximate scale factor distribution plan f of stripe pattern then _A1(x, y ₁) at coordinate (x, y ₁) locate numerical value a _RxBe a _Rx=10+0.2 * a _Max

Travel through all coordinate points of stripe pattern, try to achieve the approximate scale factor distribution plan f of stripe pattern _A1(x, y);

Step 3) is to through described step 2.1) deforming stripe figure behind the wavelet transformation removes wavelet transformation ridge error, obtains the accurate scale factor of deforming stripe image every bit, and detailed process is as follows:

Step 3.1), remove the ridge error of each coordinate points of deforming stripe image, processing procedure is:

W (x, y ₁)=W ₀+ W ₁+ W ₂+ ε ₀(5) formula

W (x, y wherein ₁) be the wavelet transformation reduced form of any coordinate points among the deforming stripe figure, wherein ε ₀Expression is led the error of bringing by the phase place second order of every bit among the deforming stripe figure, and its expression formula is:

(6) formula

Wherein, t is the intermediate computations variable,

Expression phase place second order is led, and pairing approximation scale factor differentiate is line by line obtained

Obtain

Obtain the ε of corresponding scale factor a at last ₀, obtain ε ₀Be that the one dimension plural number array of 400 row is at coordinate (x, y ₁) coordinate points on error array ε ₀(x, y ₁); W ₀, W ₁, W ₂Be respectively the simplification expression formula in the wavelet transformation computation process, form is as follows:

$W_0(a,b)=A(x,y)\frac{a\sqrt{\mathrm{\π}}}{{\left(f_b^2\mathrm{\π}\right)}^{1/4}}\exp \left({-2f}_b^2{f}_c^2{\mathrm{\π}}^2\right)$ (7) formula

(8) formula

(9) formula

Wherein,

For phase place,

Expression phase place derivative, W (x, y ₁) in pointwise deduct error ε ₀(x, y ₁), obtain the value W of accurate wavelet transformation ridge _ε(x, y ₁)=W ₀+ W ₁+ W ₂, obtain the accurate wavelet transformation ridge matrix W of each row _{400 * r}, matrix W _{400 * r}It is the complex matrix of 400 row r row;

Step 3.2), ask for the exact scale factor distribution plan f of stripe pattern _A2(x, y):

Obtain W _{400 * r}The modular matrix C of correspondence _{400 * r}(a ₁, b), search modular matrix C _{400 * r}(a ₁, the b) element of x column mean maximum, and with modular matrix C _{400 * r}(a ₁, the line label assignment of x column mean greatest member b) is given a _Amax, the approximate scale factor distribution plan f of stripe pattern then _A2(x, y ₁) at coordinate (x, y ₁) locate numerical value a _AcrFor: a _Acr=10+0.2 * a _Amax

Travel through all coordinate points of stripe pattern, try to achieve the exact scale factor distribution plan f of stripe pattern _A2(x, y);

Step 4) is done the window Fourier transform line by line to deforming stripe figure, obtains deforming stripe figure relative phase distribution plan, and detailed process is as follows:

Step 4.1), y is considered as constant, adopt the one dimension window fourier transform to deforming stripe image g (x, every row y) is handled, one dimension window fourier transform process is:

$\mathrm{WF}(b,\mathrm{\ξ})=\underset{-\∞}{\overset{+\∞}{\∫}}g\left(x\right)W_{\mathrm{\δ}}(x-b)\exp (-\mathrm{j\ξx})\mathrm{dx}$ (10) formula

Its frequency-domain expression is:

$\mathrm{WF}(f_s,y)=\underset{n=0}{\overset{n=+\&infin;}{\mathrm{\&Sigma;}}}{P}_n(f_s-{\mathrm{<figure-callout\; id="0"\; label="nf"\; filenames="HDA00003216062500031.png"\; state="\{\{state\}\}">nf</figure-callout>}}_0,y)$ (11) formula

Wherein, g (x) is delegation's stripe pattern; (ξ represents the frequency-domain calculations factor to WF for b, ξ) expression one dimension window fourier transform, and δ represents the window size factor of one dimension window Fourier, and the δ value is one dimension position of window (x-b, y) Dui Ying exact scale factor distribution plan f _A2(x, the y) value of the point on the relevant position, W _δ(x-b) expression window function, expression formula is:

$W_{\mathrm{\&delta;}}(x-b)=\frac{{\mathrm{\&delta;}}^{-1}}{\sqrt{2\mathrm{\&pi;}}}\exp [-\frac{{(b-x)}^2}{{2\mathrm{\&delta;}}^2}]$ (12) formula

Wherein, n represents order, and value is 0,1,2 successively ..., infinite; P _n(f _s-nf ₀, y) corresponding n rank frequency spectrum after the expression any point Fourier transform; f _sThe variable of expression frequency domain;

Step 4.2), to the spectral filtering after the Fourier transform and extract phase information, detailed process is as follows:

To P _n(f _s-nf ₀, y) carry out filtering and extract first order spectrum P ₁(f _s-f ₀, y), again to P ₁(f _s-f ₀, y) carry out inverse Fourier transform, obtain containing B (x, y) the exp[i φ (x of phase information, y)], calculate B (x, y) exp[i φ (x, y)] angle value can obtain containing object height information deforming stripe figure relative phase value φ (x, y), the phase value that obtains is between 0-2 π;

Travel through all coordinate points of stripe pattern, obtain the relative phase distribution plan φ of deforming stripe figure _A(x, y);

Step 5) is set up stripe pattern quality figure, adopts the flood algorithm to launch relative phase distribution plan φ _A(x y), obtains actual phase figure

Detailed process is as follows:

Step 5.1), utilize the phase gradient in the relative phase distribution plan to set up quality figure, quality figure can calculate according to following formula:

$\mathrm{\&Delta;}(x,y)=\sqrt{W_{\mathrm{rap}}^2\{{\mathrm{\&phi;}}_A(x+1,y)-{\mathrm{\&phi;}}_A(x,y)\}+W_{\mathrm{rap}}^2\{{\mathrm{\&phi;}}_A(x,y+1)-{\mathrm{\&phi;}}_A(x,y)\}}$ (13) formula

W wherein _Rap{ } is the parcel function, on dutyly it deducted or add during greater than 2 π or less than-2 π 2 π;

Step 5.2), central authorities find the highest coordinate points of mass value at the relative phase distribution plan, select this as the starting point of phase unwrapping, this point are put into the stack of a sky;

Step 5.3), judge whether stack is empty; If then the phase unwrapping process finishes, and obtains launching phase result

And enter step 6); If not, eject the point of stack top, launch not have in neighbours' point of this some the coordinate points of processing, and these untreated points are stacked;

Step 5.4), will have a few ordering in the stack according to the mass value among the quality figure, the point that quality is the highest is placed on stack top, forwards step 5.3 to) continue to handle;

Step 6) reads final expansion phase result According to classical optical grating projection from launching phase result

(x, conversion formula y) are finally tried to achieve the three-dimensional information of measuring object to object height h; Wherein conversion formula h (x y) is:

(14) formula

Wherein, l, d are the geometric parameters of measuring system, l be projector to the distance of measurement plane, d is that the CCD camera is to the distance of projector;

The expression phase changing capacity,

Be the expansion phase result,

Be initial phase result, ω ₀Angular frequency for projection grating.

Beneficial effect: compared with prior art, the present invention has the following advantages: at first, the present invention is than the three-dimensional measurement method based on wavelet transformation comparatively commonly used, and the algorithm that utilizes wavelet transformation improves the speed of asking for window size to the singularity of frequency sensitive and reduces the violent regional window of saltus step and ask for error.Lead the ridge of wavelet transformation is asked for error frequently by asking for and remove the phase place second order, can solve in the violent zone of saltus step and separate phase error problem bigger than normal, can overcome object and have the solution phase error that big saltus step and relatively gentle slope bring, obtain accurate more phase place, and improve window Fourier's noise immunity, reached measuring accuracy and stronger robustness preferably; Than Fourier's mensuration, the problem of can't local phase calculating of having avoided that global change brings makes that the phase information of calculating is more accurate; Than phase-shift method, the present invention only needs projection one secondary black and white modulation stripe image, can realize kinetic measurement; Than based on this Er Weituoke conversion improvement window Fourier three-dimensional measurement method because the multiresolution characteristic of wavelet transformation can be handled signal decomposition, can improve the accuracy rate of phase place.Secondly, the present invention obtains comparatively accurate phase information in relative phase is asked for, and makes and uses common overall flood solution pack just can obtain actual phase information accurately.At last, the present invention shows stronger robustness in processing procedure, and the ability of good actual application is arranged.To sum up, the present invention can obtain the three-dimensional elevation information of testee exactly, has robustness preferably.

Description of drawings

Fig. 1 is process flow diagram of the present invention;

Fig. 2 adopts the flood algorithm to launch the detailed process process flow diagram of phase place in the step 5;

What Fig. 3 was that CCD collects is the deforming stripe image of testee with the motorcycle backplate;

Fig. 4 is the curve map of the value delivery value behind certain some wavelet transformation among the deforming stripe figure, and horizontal ordinate represents the inverse of scale factor, and ordinate represents the mould value behind the wavelet transformation;

Fig. 5 is the mould value curve map that the value behind certain some wavelet transformation is removed the ridge error among the deforming stripe figure, and horizontal ordinate represents the inverse of scale factor, and ordinate represents the mould value of going behind the wavelet transformation after the error;

Fig. 6 is the accurate window size distribution plan that obtains;

Fig. 7 is the PHASE DISTRIBUTION figure that extracts behind the window fourier transform;

Fig. 8 utilizes the pack of flood solution to launch the absolute phase distribution plan that relative phase obtains;

Fig. 9 is the true altitude point cloud chart that obtains by PHASE DISTRIBUTION figure.

Embodiment

Below in conjunction with accompanying drawing the present invention is done further explanation.

Separate the coarse problem of phase place at the three-dimensional measurement method based on wavelet transformation in the violent zone of height saltus step, the present invention adopts Wavelet Transform to ask for the window size factor, and utilize its susceptibility to frequency factor to remove little wave crest and ask for error, improve the violent regional window size of saltus step and ask for the problem of precision, separate phase accuracy thereby improve these zones.On this basis, utilize wavelet transformation can use the characteristic of Fast Fourier Transform (FFT), shared time when improving the calculation window size, make whole process speed accelerate greatly.The present invention has not only solved the violent regional phase place of saltus step and has asked for the problem of precision, and has improved treatment effeciency.After obtaining accurate absolute phase and distributing, to the height conversion formula, finally obtain measuring the three-dimensional information of object according to the phase place of classical optical grating projection.

Under windows operating system, use VC++6.0 as programming tool, the black and white deforming stripe image that CCD collects is handled.This example adopts the plastics backplate as testee, finally obtains the more accurate whole audience absolute phase that contains the backplate three-dimensional information and distributes.

Be illustrated in figure 1 as the process flow diagram of the inventive method, the specific implementation step is as follows:

Step 1) to the testee surface, uses the CCD camera that the testee surface is taken the black and white strip image projection, obtain a panel height degree and be c, width and be r deforming stripe image g (x, y), c is that 1020, r is 1020 in the present embodiment:

G (x, y)=A (x, y)+B (x, y) cos [2 π f ₀X+ φ (x, y)] (1) formula

Wherein, (x y) is the background light distribution to A, and (x y) is the body surface reflectivity, f to B ₀Be the sine streak frequency, f in the present embodiment ₀Be that (x y) is relative phase distribution plan to be asked to 0.05, φ, (x, y) two-dimensional coordinate of expression deforming stripe image; Fig. 3 is the deforming stripe image of backplate.

Step 2), to deforming stripe figure g (x y) does wavelet transformation line by line, obtains the approximate scale factor of deforming stripe image every bit, and detailed process is as follows:

Step 2.1), y is considered as constant y ₁, adopt one-dimensional wavelet transform to deforming stripe image g (x, y y) ₁Row is handled, and processing procedure is:

${W_y}_1(a_1,b)={\&Integral;}_{-\&infin;}^{+\&infin;}g(x,y_1){M^{*}}_{a,b}\left(x\right)\mathrm{dx}$ (2) formula

Wherein, y ₁The row number of representing certain delegation; A is scale factor, and span is 10 to 90, gets a value every 0.2; B is shift factor, and span is 1 to the width r of striped, gets a value every 1, and unit is pixel; Obtain through one-dimensional wavelet transform

Be the two-dimentional complex matrix of one 400 capable r row, claim y ₁The row wavelet transform matrix, a ₁It is matrix

The line label of middle element; Wherein,

M ^* _{A, b}(x) be M _{A, b}(x) conjugate function; Wavelet function M (x) according to formula (3) obtains M ^* _{A, b}(x) expression formula is suc as formula (4):

$M\left(x\right)=\frac{1}{{\left(f_b^2\mathrm{\&pi;}\right)}^{1/4}}\exp \left(2\mathrm{\&pi;}{\mathrm{if}}_{c}x\right)\exp \left(\frac{{-x}^2}{{2f}_b^2}\right)$ (3) formula

${M^{*}}_{a,b}\left(x\right)=\frac{1}{\sqrt{a}{\left(f_b^2\mathrm{\&pi;}\right)}^{1/4}}\exp \left({-2\mathrm{\&pi;if}}_c\frac{x-b}{a}\right)\exp [-\frac{{(x-b)}^2}{{2f}_b^2{a}^2}];$ (4) formula

Wherein, f _bBe the bandwidth of wavelet function, value is 0.8, f in the present embodiment _cBe the centre frequency of wavelet function, value is that 1, i is complex unit in the present embodiment, and x is independent variable.

Step 2.2), ask for stripe pattern g (x, approximate scale factor distribution plan f y) _A1(x, y ₁):

Obtain two-dimentional complex matrix The modular matrix of correspondence

Its form is:

(5) formula

The a single point experiment effect is searched for modular matrix as shown in Figure 4

The element of x column mean maximum, and with modular matrix The line label assignment of x column mean greatest member give a _Max, the approximate scale factor distribution plan f of stripe pattern then _A1(x, y ₁) at coordinate (x, y ₁) locate numerical value a _RxBe a _Rx=10+0.2 * a _MaxTravel through all coordinate points of stripe pattern, try to achieve the approximate scale factor distribution plan f of stripe pattern _A1(x, y) as shown in Figure 4.

Step 3) is to through step 2.1) deforming stripe figure behind the wavelet transformation removes wavelet transformation ridge error, obtains the accurate scale factor of deforming stripe image every bit, and detailed process is as follows:

Step 3.1), remove through step 3.1) the ridge error of each coordinate points of deforming stripe image behind the wavelet transformation, because small echo ridge method is only chosen

Value before, and be easy to occur the violent zone of phase hit in actual use,

Also no longer be to go to zero and the amount that can eliminate, exist usually

Change big zone To tend towards stability, at this moment

To ridge to ask for precision influence very low, so can suppose in the violent zone of saltus step

Be 0.Therefore we can be shown phase meter:

(6) formula

Even in processing procedure, occur violent hop region like this, owing to considered The error of bringing will significantly reduce the error of calculation, so the simplification process is in the wavelet transform process process:

W (x, y ₁)=W ₀+ W ₁+ W ₂+ ε ₀(7) formula

W ₀, W ₁, W ₂Be respectively the simplification expression formula in the wavelet transformation computation process, form is as follows:

$W_0(a,b)=A(x,y)\frac{a\sqrt{\mathrm{\&pi;}}}{{\left(f_b^2\mathrm{\&pi;}\right)}^{1/4}}\exp (-2f_b^2{f}_c^2{\mathrm{\&pi;}}^2)$ (8) formula

(9) formula

(10) formula

W (x, y wherein ₁) be the wavelet transformation reduced form of any pixel among the deforming stripe figure, wherein ε ₀Expression is owing to the phase place second order of every bit among the deforming stripe figure is led the error of bringing, and its expression formula is:

(11) formula

Wherein t is the intermediate computations variable, Expression phase place second order is led, Expression phase place single order is led, and wherein comprises variable a, as seen

When big,

To depart from ridge frequently.Because wavelet transformation is constantly to use different scale factor a to ask for ridge frequently (wherein scale factor a changes to 90 with 0.2 amplitude from 10) to every bit, so wavelet transformation W (x, y ₁)=W ₀+ W ₁+ W ₂+ ε ₀In each point by frequently removing ε ₀, ask for W again _ε(x, y ₁)=W ₀+ W ₁+ W ₂Ridge, just can obtain comparatively accurate

ε ₀In have only

A unknown quantity utilizes following method to be similar to and asks for

Consider

After, f exist one with

Relevant deviation is made as

(12) formula

Differentiate gets to f

(13) formula

Pairing approximation scale factor differentiate is line by line obtained

Obtain Obtain the ε of corresponding scale factor a at last ₀, obtain ε ₀Be that the one dimension plural number array of 400 row is at coordinate (x, y ₁) error array ε ₀(x, y ₁).W (x, y ₁) in pointwise deduct error ε ₀(x, y ₁), obtain the value W of accurate wavelet transformation ridge _ε(x, y ₁)=W ₀+ W ₁+ W ₂, obtain the accurate wavelet transformation ridge matrix W of each row _{400 * r}, matrix W _{400 * r}It is the complex matrix of 400 row r row.

Step 3.2), ask for the exact scale factor distribution plan f of stripe pattern _A2(x, y):

Obtain W _{400 * r}The modular matrix C of correspondence _{400 * r}(a ₁, b), search modular matrix C ₄₀₀* r (a ₁, the b) element of x column mean maximum, and with modular matrix C _{400 * r}(a ₁, the line label assignment of x column mean greatest member b) is given a _Amax, the approximate scale factor distribution plan f of stripe pattern then _A2(x, y ₁) at coordinate (x, y ₁) locate numerical value a _AcrFor: a _Acr=10+0.2 * a _Amax

Travel through all coordinate points of stripe pattern, try to achieve the exact scale factor distribution plan f of stripe pattern _A2(x, y) as shown in Figure 6.

Step 4) is done the window Fourier transform line by line to deforming stripe figure, obtains deforming stripe figure relative phase distribution plan, and detailed process is as follows:

Step 4.1), y is considered as constant, adopt the one dimension window fourier transform to deforming stripe image g (x, every row y) is handled, one dimension window fourier transform process is:

$\mathrm{WF}(b,\mathrm{\&xi;})=\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}g\left(x\right)W_{\mathrm{\&delta;}}(x-b)\exp (-\mathrm{j\&xi;x})\mathrm{dx}$ (14) formula

Its frequency-domain expression is:

$\mathrm{WF}(f_s,y)=\underset{n=0}{\overset{n=+\&infin;}{\mathrm{\&Sigma;}}}{P}_n(f_s-{\mathrm{<figure-callout\; id="0"\; label="nf"\; filenames="HDA00003216062500031.png"\; state="\{\{state\}\}">nf</figure-callout>}}_0,y)$ (15) formula

Wherein, g (x) is delegation's stripe pattern; (ξ represents the frequency-domain calculations factor to WF for b, ξ) expression one dimension window fourier transform, and δ represents the window size factor of one dimension window Fourier, and the δ value is one dimension position of window (x-b, y) Dui Ying exact scale factor distribution plan f _A2(x, the y) value of the point on the relevant position, W _δ(x-b) expression window function, expression formula is:

$W_{\mathrm{\&delta;}}(x-b)=\frac{{\mathrm{\&delta;}}^{-1}}{\sqrt{2\mathrm{\&pi;}}}\exp [-\frac{{(b-x)}^2}{{2\mathrm{\&delta;}}^2}]$ (16) formula

Wherein, n represents order, and value is 0,1,2 successively ..., infinite; P _n(f _s-nf ₀, y) corresponding n rank frequency spectrum after the expression any point Fourier transform; f _sThe variable of expression frequency domain;

Step 4.2), to the spectral filtering after the Fourier transform and extract phase information, detailed process is as follows:

To P _n(f _s-nf ₀, y) carry out filtering and extract first order spectrum P ₁(f _s-f ₀, y), again to P ₁(f _s-f ₀, y) carry out inverse Fourier transform, obtain containing B (x, y) the exp[i φ (x of phase information, y)], calculate B (x, y) exp[i φ (x, y)] angle value can obtain containing object height information deforming stripe figure relative phase value φ (x, y), the phase value that obtains is between 0-2 π;

Travel through all coordinate points of stripe pattern, obtain the relative phase distribution plan φ of deforming stripe figure _A(x, y), the wrapped phase lab diagram as shown in Figure 7.

Step 5) is set up stripe pattern quality figure, adopts the flood algorithm to launch relative phase distribution plan φ _A(x y), obtains actual phase figure

Detailed process is as follows:

Step 5.1), utilize the phase gradient among the relative phase figure to set up quality figure, the reliability that phase gradient also can indicate this place's phase place to find the solution, phase gradient is more big, illustrate that this zone may exist object height discontinuous or find the solution problems such as mistake, reliability is lower, and phase gradient is more little, illustrate that this zone should belong to mild zone, reliability is higher; Quality figure can calculate according to following formula:

$\mathrm{\&Delta;}(x,y)=\sqrt{W_{\mathrm{rap}}^2\{{\mathrm{\&phi;}}_A(x+1,y)-{\mathrm{\&phi;}}_A(x,y)\}+W_{\mathrm{rap}}^2\{{\mathrm{\&phi;}}_A(x,y+1)-{\mathrm{\&phi;}}_A(x,y)\}}$ (17) formula

W wherein _Rap{ } is the parcel function, on dutyly it deducted or add during greater than 2 π or less than-2 π 2 π.

Step 5.2), central authorities find the highest coordinate points of mass value at the relative phase distribution plan, select this as the starting point of phase unwrapping, this point are put into the stack of a sky;

Step 5.3), judge whether stack is empty; If then the phase unwrapping process finishes, and obtains launching phase result And enter step 6); If not, eject the point of stack top, launch not have in neighbours' point of this some the coordinate points of processing, and these untreated points are stacked.

Step 5.4), will have a few ordering in the stack according to the mass value among the quality figure, the point that quality is the highest is placed on stack top, forwards step 5.3 to) continue to handle, final phase unwrapping lab diagram is as shown in Figure 8.

Step 6) reads final expansion phase result

According to classical optical grating projection from launching phase result

(x, conversion formula y) are finally tried to achieve the three-dimensional information of measuring object to object height h; Wherein conversion formula h (x y) is:

(18) formula

Wherein, l, d are the geometric parameters of measuring system, l be projector to the distance of measurement plane, l is 164 centimetres in the present embodiment, d be the CCD camera to the distance of projector, d is 58 centimetres in the present embodiment;

The expression phase changing capacity,

Be the expansion phase result,

Be initial phase result, ω ₀Angular frequency for projection grating.Fig. 9 is for measuring object dimensional information.

The above only is preferred implementation of the present invention; should be pointed out that for those skilled in the art, under the prerequisite that does not break away from the principle of the invention; can also make some improvements and modifications, these improvements and modifications also should be considered as protection scope of the present invention.

Claims (1)

1. the improvement window Fourier three-dimensional measurement method based on wavelet transformation is characterized in that, comprises the steps:

Step 1) to the testee surface, uses the CCD camera that the testee surface is taken the black and white strip image projection, obtain a panel height degree and be c, width and be r deforming stripe image g (x, y):

G (x, y)=A (x, y)+B (x, y) cos[2 π f ₀X+ φ (x, y)] (1) formula

Wherein, (x y) is the background light distribution to A, and (x y) is the body surface reflectivity, f to B ₀Be the sine streak frequency, (x y) is relative phase distribution plan to be asked to φ, (x, y) two-dimensional coordinate of expression deforming stripe image;

Step 2), to described deforming stripe figure g (x y) does wavelet transformation line by line, obtains the approximate scale factor of deforming stripe image every bit, and detailed process is as follows:

Step 2.1), y is considered as constant y ₁, adopt one-dimensional wavelet transform to described deforming stripe image g (x, y y) ₁Row is handled, and processing procedure is:

${W_y}_1(a_1,b)={\&Integral;}_{-\&infin;}^{+\&infin;}g(x,y_1){M^{*}}_{a,b}\left(x\right)\mathrm{dx}$ (2) formula

Wherein, y ₁The row number of representing certain delegation; A is scale factor, and span is 10 to 90, gets a value every 0.2; B is shift factor, and span is 1 to the width r of striped, gets a value every 1, and unit is pixel; Obtain through one-dimensional wavelet transform

Be the two-dimentional complex matrix of one 400 capable r row, claim y ₁The row wavelet transform matrix, a ₁It is matrix

The line label of middle element; Wherein,

M ^* _{A, b}(x) be M _{A, b}(x) conjugate function; Wavelet function M (x) according to formula (3) obtains M ^* _{A, b}(x) expression formula is suc as formula (4):

$M\left(x\right)=\frac{1}{{\left(f_b^2\mathrm{\&pi;}\right)}^{1/4}}\exp \left(2{\mathrm{\&pi;if}}_{c}x\right)\exp \left(\frac{{-x}^2}{{2f}_b^2}\right)$ (3) formula

${M^{*}}_{a,b}\left(x\right)=\frac{1}{\sqrt{a}{\left(f_b^2\mathrm{\&pi;}\right)}^{1/4}}\exp (-{2\mathrm{\&pi;if}}_c\frac{x-b}{a})\exp [-\frac{{(x-b)}^2}{{2f}_b^2{a}^2}];$ (4) formula

Wherein, f _bBe the bandwidth of wavelet function, f _cBe the centre frequency of wavelet function, i is complex unit, and x is independent variable;

Step 2.2), ask for stripe pattern g (x, approximate scale factor distribution plan f y) _A1(x, y ₁):

Calculate two-dimentional complex matrix

The modular matrix of correspondence

The search modular matrix

The element of x column mean maximum, and with modular matrix

The line label assignment of x column mean greatest member give a _Max, the approximate scale factor distribution plan f of stripe pattern then _A1(x, y ₁) at coordinate (x, y ₁) locate numerical value a _RxBe a _Rx=10+0.2 * a _Max

Travel through all coordinate points of stripe pattern, try to achieve the approximate scale factor distribution plan f of stripe pattern _A1(x, y);

Step 3) is to through described step 2.1) deforming stripe figure behind the wavelet transformation removes wavelet transformation ridge error, obtains the accurate scale factor of deforming stripe image every bit, and detailed process is as follows:

Step 3.1), remove the ridge error of each coordinate points of deforming stripe image, processing procedure is:

W (x, y ₁)=W ₀+ W ₁+ W ₂+ ε ₀(5) formula

W (x, y wherein ₁) be the wavelet transformation reduced form of any coordinate points among the deforming stripe figure, wherein ε ₀Expression is led the error of bringing by the phase place second order of every bit among the deforming stripe figure, and its expression formula is:

(6) formula

Wherein, t is the intermediate computations variable,

Expression phase place second order is led, and pairing approximation scale factor differentiate is line by line obtained

Obtain

Obtain the ε of corresponding scale factor a at last ₀, obtain ε ₀Be that the one dimension plural number array of 400 row is at coordinate (x, y ₁) coordinate points on error array ε ₀(x, y ₁); W ₀, W ₁, W ₂Be respectively the simplification expression formula in the wavelet transformation computation process, form is as follows:

$W_0(a,b)=A(x,y)\frac{a\sqrt{\mathrm{\&pi;}}}{{\left(f_b^2\mathrm{\&pi;}\right)}^{1/4}}\exp \left({-2f}_b^2{f}_c^2{\mathrm{\&pi;}}^2\right)$ (7) formula

(8) formula

(9) formula

Wherein,

For phase place,

Expression phase place derivative, W (x, y ₁) in pointwise deduct error ε ₀(x, y ₁), obtain the value W of accurate wavelet transformation ridge _ε(x, y ₁)=W ₀+ W ₁+ W ₂, obtain the accurate wavelet transformation ridge matrix W of each row _{400 * r}, matrix W _{400 * r}It is the complex matrix of 400 row r row;

Step 3.2), ask for the exact scale factor distribution plan f of stripe pattern _A2(x, y):

Obtain W _{400 * r}The modular matrix C of correspondence _{400 * r}(a ₁, b), search modular matrix C _{400 * r}(a ₁, the b) element of x column mean maximum, and with modular matrix C _{400 * r}(a ₁, the line label assignment of x column mean greatest member b) is given a _Amax, the approximate scale factor distribution plan f of stripe pattern then _A2(x, y ₁) at coordinate (x, y ₁) locate numerical value a _AcrFor: a _Acr=10+0.2 * a _Amax

Travel through all coordinate points of stripe pattern, try to achieve the exact scale factor distribution plan f of stripe pattern _A2(x, y);

Step 4) is done the window Fourier transform line by line to deforming stripe figure, obtains deforming stripe figure relative phase distribution plan, and detailed process is as follows:

Step 4.1), y is considered as constant, adopt the one dimension window fourier transform to deforming stripe image g (x, every row y) is handled, one dimension window fourier transform process is:

$\mathrm{WF}(b,\mathrm{\&xi;})=\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}g\left(x\right)W_{\mathrm{\&delta;}}(x-b)\exp (-\mathrm{j\&xi;x})\mathrm{dx}$ (10) formula

Its frequency-domain expression is:

$\mathrm{WF}(f_s,y)=\underset{n=0}{\overset{n=+\&infin;}{\mathrm{\&Sigma;}}}{P}_n(f_s-{\mathrm{nf}}_0,y)$ (11) formula

Wherein, g (x) is delegation's stripe pattern; (ξ represents the frequency-domain calculations factor to WF for b, ξ) expression one dimension window fourier transform, and δ represents the window size factor of one dimension window Fourier, and the δ value is one dimension position of window (x-b, y) Dui Ying exact scale factor distribution plan f _A2(x, the y) value of the point on the relevant position, W _δ(x-b) expression window function, expression formula is:

$W_{\mathrm{\&delta;}}(x-b)=\frac{{\mathrm{\&delta;}}^{-1}}{\sqrt{2\mathrm{\&pi;}}}\exp [-\frac{{(b-x)}^2}{{2\mathrm{\&delta;}}^2}]$ (12) formula

Wherein, n represents order, and value is 0,1,2 successively ..., infinite; P _n(f _s-nf ₀, y) corresponding n rank frequency spectrum after the expression any point Fourier transform; f _sThe variable of expression frequency domain;

Step 4.2), to the spectral filtering after the Fourier transform and extract phase information, detailed process is as follows:

To P _n(f _s-nf ₀, y) carry out filtering and extract first order spectrum P ₁(f _s-f ₀, y), again to P ₁(f _s-f ₀, y) carry out inverse Fourier transform, obtain containing B (x, y) the exp[i φ (x of phase information, y)], calculate B (x, y) exp[i φ (x, y)] angle value can obtain containing object height information deforming stripe figure relative phase value φ (x, y), the phase value that obtains is between 0-2 π;

Travel through all coordinate points of stripe pattern, obtain the relative phase distribution plan φ of deforming stripe figure _A(x, y);

Step 5) is set up stripe pattern quality figure, adopts the flood algorithm to launch relative phase distribution plan φ _A(x y), obtains actual phase figure

Detailed process is as follows:

Step 5.1), utilize the phase gradient in the relative phase distribution plan to set up quality figure, quality figure can calculate according to following formula:

$\mathrm{\&Delta;}(x,y)=\sqrt{W_{\mathrm{rap}}^2\{{\mathrm{\&phi;}}_A(x+1,y)-{\mathrm{\&phi;}}_A(x,y)\}+W_{\mathrm{rap}}^2\{{\mathrm{\&phi;}}_A(x,y+1)-{\mathrm{\&phi;}}_A(x,y)\}}$ (13) formula

W wherein _Rap{ } is the parcel function, on dutyly it deducted or add during greater than 2 π or less than-2 π 2 π;

Step 5.2), central authorities find the highest coordinate points of mass value at the relative phase distribution plan, select this as the starting point of phase unwrapping, this point are put into the stack of a sky;

Step 5.3), judge whether stack is empty; If then the phase unwrapping process finishes, and obtains launching phase result And enter step 6); If not, eject the point of stack top, launch not have in neighbours' point of this some the coordinate points of processing, and these untreated points are stacked;

Step 5.4), will have a few ordering in the stack according to the mass value among the quality figure, the point that quality is the highest is placed on stack top, forwards step 5.3 to) continue to handle;

Step 6) reads final expansion phase result

According to classical optical grating projection from launching phase result

(x, conversion formula y) are finally tried to achieve the three-dimensional information of measuring object to object height h; Wherein conversion formula h (x y) is:

(14) formula

Wherein, l, d are the geometric parameters of measuring system, l be projector to the distance of measurement plane, d is that the CCD camera is to the distance of projector;

The expression phase changing capacity,

Be the expansion phase result,

Be initial phase result, ω ₀Angular frequency for projection grating.

Images