CN108680119A - A kind of subregional single width fast phase method of deploying - Google Patents
A kind of subregional single width fast phase method of deploying Download PDFInfo
- Publication number
- CN108680119A CN108680119A CN201810352746.7A CN201810352746A CN108680119A CN 108680119 A CN108680119 A CN 108680119A CN 201810352746 A CN201810352746 A CN 201810352746A CN 108680119 A CN108680119 A CN 108680119A
- Authority
- CN
- China
- Prior art keywords
- region
- phase
- subregional
- result
- single width
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B11/00—Measuring arrangements characterised by the use of optical techniques
- G01B11/24—Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures
- G01B11/25—Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures by projecting a pattern, e.g. one or more lines, moiré fringes on the object
- G01B11/2518—Projection by scanning of the object
- G01B11/2527—Projection by scanning of the object with phase change by in-plane movement of the patern
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B11/00—Measuring arrangements characterised by the use of optical techniques
- G01B11/24—Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures
- G01B11/25—Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures by projecting a pattern, e.g. one or more lines, moiré fringes on the object
- G01B11/254—Projection of a pattern, viewing through a pattern, e.g. moiré
Landscapes
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Image Analysis (AREA)
- Image Processing (AREA)
Abstract
The invention discloses a kind of subregional single width fast phase method of deploying.Self-adapting multi-dimension decomposition is carried out to bar graph first with the Empirical mode decomposition (BSEMD) of two-dimentional sinusoidal auxiliary, obtained high frequency result component is for detecting object edge in striped, obtained middle stripe carrier component obtains wrapped phase and instantaneous frequency by Hilbert transform, and obtained low frequency component is identified for bar graph shadow region；Then edge and shadow region are merged and carries out morphological operation and form low quality region；High quality region is unfolded followed by the mode that flood is spread, the low quality region phase unwrapping of Quality Map guidance extension is established with instantaneous frequency；Finally each region phase is merged to form complete continuous phase.For the present invention by dividing rapid deployment of the region using different phase unwrapping strategy realization wrapped phases, entire method is efficient and precision is high, has larger application prospect in single width striped three-dimensional Fast measurement system.
Description
Technical field
The invention belongs to three dimension reconstruction technical fields, more particularly to a kind of subregional single width fast phase expansion
Method.
Background technology
The advantages that fringe projection measuring technique is non-contact with its, quick, accurate has become a kind of three-dimensional being very welcomed
Surface shape measurement technology is all widely used in national defense safety, advanced manufacture, reverse-engineering and historical relic's protection field.Wherein, more
The optical measuring technique of spoke line can realize the measurement of high-precision and less constraint.However, this method needs to keep
The static projection multiple image of object cannot achieve dynamic 3 D measurement, therefore single width fringe projection technology is come into being.
Fringe projection technology shoots the striped modulated by body surface by projecting sinusoidal grating striped to testee
The three-dimensional information of pattern, body surface is included in fringe phase, it is therefore desirable to be demodulated to stripe pattern.Based on single width item
The phase acquiring technology of line projection mainly has：Fourier transformation, window Fourier transform, S-transformation, wavelet transformation and two-dimensional empirical
Mode decomposition (BSEMD) etc..Wherein, BSEMD is as a kind of new adaptive Time Frequency Analysis method, in analysis unstable signal
When, due to being driven by own signal, possess has better time frequency analysis ability than other methods, and application scenarios of the invention are just
It is based on BSEMD.
These phase demodulatings based on transform domain are obtained by contact transformation anyway, to obtain at main value phase
In range [- π, π), it is therefore desirable to be launched into continuous phase.Phase unwrapping only needs more adjacent 2 points packet under normal circumstances
Wrap up in whether phase difference is completed more than π to compensate 2 π offsets.However, due to lack sampling, noise, invalid shadow region, object
The bad phase of formation, tradition is sharply waited to be easy to form " bracing wire " phenomenon by row method of deploying in surface.
Quality Map method is the most common type in phase developing method, and Quality Map describes each pixel phase in bar graph
Quality is to instruct phase integral to be unfolded.The path method of deploying of Quality Map guiding, although the propagation of phase error can be controlled,
It is that it compares the quality of each point to select high quality point to be preferentially unfolded, it is long that this results in Riming time of algorithm, greatly
The big rapid charater for reducing single width fringe projection.
Invention content
In order to solve the technical issues of above-mentioned background technology proposes, the present invention is intended to provide a kind of subregional single width is quick
Phase developing method is unfolded phase by subregion, reduces the run time of Quality Map method.
In order to achieve the above technical purposes, the technical scheme is that：
A kind of subregional single width fast phase method of deploying, includes the following steps：
(1) it uses BSEMD to carry out multi-resolution decomposition to bar graph I (x, y), obtains the exploded view BIMF of k width different scalesi
(x,y)：
Wherein, BIMFiFor intrinsic mode function, r (x, y) is surplus, and determines parameter k1、k2To divide high fdrequency componentThe striped carrier component of enhancingAnd low frequency component
(2) Hilbert transform is carried out to striped carrier component and obtains wrapped phase and instantaneous frequency fins(x, y), and profit
With instantaneous frequency fins(x, y) establishes Quality Map Q (x, y)；
(3) edge detection is carried out to high fdrequency component, shadow region identification is carried out to low frequency component, both fusions are as a result, again
Result applied morphology after fusion is operated to form low quality region, remaining region is high quality region；
(4) label distinguishes each connected domain, using flood sprawling expansion phase in each high quality region, extension it is low
Phase is unfolded using Quality Map method in quality region；
(5) adjacent high and low quality region phase merges, until ultimately forming complete continuous phase.
Further, in step (2), Quality Map Q (x, y) is obtained according to the following formula：
Q (x, y)=| fins(x,y)-f0|
Wherein, f0For striped carrier frequency.
Further, in step (3), edge detection is carried out to high fdrequency component using Canny operators, by two values matrix
Edge pixel values are set to 1 and are used as edge detection results in MASK1.
Further, in step (3), shadow region is carried out to low frequency component and knows method for distinguishing：Low frequency component is used
Otsu methods carry out multi-threshold segmentation, and as a result the minimum pixel forming region of middle gray value is shadow region, by two values matrix MASK2
Middle shadow region pixel value is set to 1 and is used as shadow region recognition result.
Further, in step (3), the fusion refers to obtaining two values matrix MASK1 and MASK2 progress or operation
To MASK3.
Further, in step (3), the morphological operation is as follows：
Wherein, BI is circular structural element, the pixel groups that the low quality region is 1 by two values matrix MASK intermediate values
At.
Further, in step (4), the step of marking each connected domain, is as follows：
(401) each 4 in the region that the quick recognition matrix MASK values of MATLAB tool functions bwlabel are 0 are used to connect
Lead to region, grey scale pixel value is respectively 1,2 in each region ..., and m, wherein m are high quality areal；
(402) each 4 in the region that the quick recognition matrix MASK values of MATLAB tool functions bwlabel are 1 are used to connect
Lead to region, grey scale pixel value is that respectively m+1, m+2 ..., m+n, wherein n are low quality areal in each region；
(403) result of step (401) and (402) is done into matrix and operation, operating result is each connected component labeling
As a result, being indicated with matrix MARK, MARK values are 1,2 ..., and m indicates that each high quality region, MARK values are m+1, m+ respectively
2 ..., m+n indicates each low quality region respectively.
Further, in step (4), the meaning in the low quality region of the extension：Low quality region in matrix MARK
Including boundary pixel point, and boundary pixel point is uniformly selected in adjacent high quality region.
Further, in step (5), two neighboring region phase merging process is indicated with following formula：
Wherein, KjFor phase amalgamation result, round () is the operation that rounds up,It is boundary element phase in high quality
Expansion is as a result, ψ in regionjIt is unfolded as a result, n is boundary pixel sum in the low quality region of extension for boundary element.
The advantageous effect brought using above-mentioned technical proposal：
(1) while the present invention obtains wrapped phase using BSEMD, Quality Map is established using intermediate result instantaneous frequency,
Greatly reduce the time established used in Quality Map.
(2) present invention is generated using adaptive region, and region division depends on physical surface characteristics and lighting angle, compares
Possess more flexibilities in fixed region partitioning method, it is more scientific and reasonable to the processing of zone boundary.
(3) present invention is unfolded compared to global Quality Map method, possesses faster development rate so that single width striped is more
It is measured suitable for quick three-dimensional.
Description of the drawings
Fig. 1 is the overall flow figure of the present invention；
Fig. 2 is embodiment bar graph；
Fig. 3 is BSEMD algorithm flow charts；
Fig. 4 is embodiment wrapped phase figure；
Fig. 5 is the Quality Map that embodiment is established；
Fig. 6 is embodiment edge detection schematic diagram；
Fig. 7 is embodiment shadow region testing result figure；
Fig. 8 is embodiment region recognition number schematic diagram；
Fig. 9, Figure 10 are that embodiment boundary element aids in illustrating schematic diagram；
Figure 11 is the final phase unwrapping result figure of embodiment.
Specific implementation mode
Below with reference to attached drawing, technical scheme of the present invention is described in detail.
As shown in Figure 1, a kind of subregional single width fast phase method of deploying that the present invention designs, just first with two dimension
The Empirical mode decomposition (BSEMD) of string auxiliary carries out self-adapting multi-dimension decomposition, obtained high frequency result component to bar graph
For detecting object edge in striped, obtained middle stripe carrier component obtains wrapped phase and wink by Hilbert transform
When frequency, obtained low frequency component for bar graph shadow region identify.Then, edge and shadow region are merged and carries out shape
State operation forms low quality region.High quality region is unfolded followed by the mode that flood is spread, matter is established with instantaneous frequency
Spirogram instructs low quality region phase unwrapping.Finally each region phase is merged to form complete continuous phase.This is described below
The detailed process of method.
Step 1：It is decomposed by the way that BSEMD algorithms are (as shown in Figure 2) to bar graph Im (x, y), high frequency division is divided from striped
AmountStriped carrier component PM, the low frequency component of enhancingThe algorithm flow chart is as shown in figure 3, will grasp
It is described as follows as step：
1) init attributes amount：
Ii(x, y)=Im (x, y), each exploded view BIMFiThe ordinal number i=1, c, l of (x, y) are striped broadband and height, k1
=0；
BIMFiThe statistical average period AvDistPre=1, each BIMF of (x, y)iThe set Exe in (x, y) statistical average period
For empty set；
BIMFiENERGY E=0 of (x, y), each BIMFi(x, y) energy quantity set Ep is empty set；
BIMFiThe global frequencies fR=0, each BIMF of (x, y)i(x, y) global frequencies integrate fRw as empty set；
FRk is empty set, with storageEach BIMF lateriThe global frequencies of (x, y)；
ARk is empty set, with storageEach BIMF lateriThe global amplitude of (x, y)；
Rfk is empty set, with storageThe global frequencies ratio of two neighboring BIMF later；
Final striped carrier component is sig=0.
2) use BSEMD methods to Ii(x, y) carries out single layer Scale Decomposition, and the wherein input of BSEMD single layers analytic function is Ii
(x, y) and statistical average cycle parameter AvDistPre, exports as BIMFi(x, y), BSEMD algorithms are as follows：
2.1) to Ii(x, y) uses Fast Experience Mode Decomposition (Enhanced Fast Empirical Mode
Decomposition, EFEMD) carry out single layer scale decomposition, wherein the input of EFEMD single layers analytic function be Ii(x, y) with
And statistical average cycle parameter AvDistPre, it exports as BIMFout(x, y) and new statistical average period AvDistPre0；If
It setsFor IiThe initial decomposition figure of (x, y), specifically, the single layer decomposition step of EFEMD is such as
Under：
2.1.1) to Ii(x, y) carries out minute surface continuation, generates Ie(x,y)；Specifically, the minute surface continuation step based on Matlab
It is rapid as follows：
Determine that continuation size is es=round [min (c, l)/10]；
It extracts image Im (x, y) width and is the outline border of es, and ask its minute surface as peripheral continuation value, i.e.,：Horizontal top is prolonged
It is ext1=flipud [Im (1 to open up value:es,1:C)], wherein flipud () is that the matrix carried in Matlab spins upside down letter
Number；Horizontal bottom continuation value is ext2=flipud [Im (l-es+1:l,1:c)]；Left side continuation value is ext3=fliplr [Im
(1:l,1:Es)], wherein fliplr () is matrix carried in Matlab or so overturning function；The right continuation value is ext4=
fliplr[Im(1:l,c-es+1:c)]；Upper left corner continuation value is ext5=fliplr { rot90 [Im (1:es,1:Es)] },
Middle rot90 () is for Matlab from tape function to realize that matrix is rotated by 90 ° counterclockwise；Upper right corner continuation value is ext6=Im (1:
es,c-es+1:c)T, ()TFor Matrix Calculating transposition；Lower left corner continuation value is ext7=Im (l-es+1:l,1:es)T；The lower right corner
Continuation value is ext8=fliplr { rot90 [Im (l-es+1:l,c-es+1:c)]}；
All continuation values of above-mentioned steps and Im (x, y) are combined, finally obtaining the result after minute surface continuation is：
2.1.2) to Ie(x, y) carries out morphological dilation twice respectively, i.e.,：
MaxD (x, y)=imdilate (Ie(x,y),ms)
MinD (x, y)=- imdilate (- Ie(x,y),ms)
Wherein ms=strel (' disk ', AvDistPre), strel () are to be used from tape function in the tool boxes Matlab
To create the structural element of designated shape, imdilate () is the morphological dilations function carried in Matlab；Seek Ie(x,y)
Extreme maximum distribution be：
MaxM (x, y)=~(Ie(x,y)-MaxD(x,y))
IeThe minimum of (y, x) is distributed as：
MinM (x, y)=~(Ie(x,y)-MinD(x,y))
Wherein "~" is " logic NOT " operation, i.e., the value that element in matrix is 0 is set to 1, non-zero value and is set to 0；
2.1.3 the mean space distance of all extreme values) is calculated：
Wherein round () operates for round numbers；If AvDistNew>AvDistPre enables AvDistPre=
Otherwise AvDistNew enables AvDistPre=2 × AvDistPre；
2.1.4) two morphological dilations results in step 2.1.2 are averaging, and the average value to obtaining carries out two
Tie up convolution smooth operation：
Wherein mC=fspecial (" disk ", round (AvDistPre/2)), conv2 () are in the tool boxes Matlab
Two-dimensional convolution handling function, fspecial () be the tool boxes Matlab in establish predefine filter operator function.
2.1.5) by the sharpening result obtained in step 2.1.4 from IeIt is subtracted in (x, y)：
B (x, y)=Ie(x,y)-SM(x,y)
The value in continuation region in B (x, y) is removed according to step 2.1.1, obtained result BIMFout(x, y) is to use
One layer of exploded view that EFEMD is obtained, EFEMD final outputs BIMFout(x, y) and AvDistPre0, noteFor initial decomposition figure；
2.2) obtained according to step 2.1Designing two dimension Sine distribution s (x, y) is：
S (x, y)=acos (2 π fy) cos (2 π fx)
Wherein, amplitude a is designed as the maximum value of exploded view two-dimensional matrix all elements：
Max () operates for maximizing；Frequency f can be acquired by following formula：
2.3) by s (x, y) respectively withMutually adduction is subtracted each other, i.e.,：
It is rightEFEMD single layer decomposition is carried out, input parameter is still AvDistPre, obtains the first katolysis
FigureThe new parameter AvDistPre exported simultaneouslys1；Similarly, rightCarry out EFEMD single layers
(input parameter is still AvDistPre) is decomposed, the second katolysis figure is obtainedAnd newly exported
AvDistPres2；
2.4) to obtainedWithIt averages：
BIMFi(x, y) is the final exploded view of image Im (x, y)；Then to two output parameter AvDistPres1With
AvDistPres2It averages, enables：
Wherein round () operates for round numbers, and further, extension Exe is Exe=[Exe AvDistPrei], it is right
AvDistPreiIt is corrected as follows：
If (AvDistPrei＜ 2Exe (end))
Wherein Exe (end) refers to Exe and concentrates the last one value；By BIMFi(x, y) is from IiIt subtracts, enables in (x, y)：Ii+1(x,
Y)=Ii(x,y)-BIMFi(x, y), Ii+1(x, y) is as the input picture decomposed next time.
3) BIMF is calculatediThe two-dimensional autocorrelation function of (x, y)：
Wherein (τ1,τ2) be (x, y) corresponding direction when seeking auto-correlation function time delay, and 0≤τ1≤ 2x-1,0≤τ2≤
2y-1；
According to Ri(τ1,τ2), estimate BIMFiThe energy value of (x, y) is approximately：
E=max [Ri(τ1,τ2)]
Max [] operates for maximizing；Ep is extended to Ep=[Ep E]；
According to Ri(τ1,τ2), estimate BIMFiThe global frequencies fR of (x, y), steps are as follows：
3.1) R is extractedi(τ1,τ2) l rows data and take the row data to be located at the right half part of non-negative axis, i.e. Rl=Ri
(l,c:End), find out all extreme points (both ends boundary point is not counting extreme point) of Rl and find first extreme point and be located at Rl's
Coordinate value Rld, Rld are approximately BIMFiThe period of (x, y) main component in the row direction is arranged if Rl does not have extreme point
Rld=1；
Similarly, R is extractedi(τ1,τ2) c column datas and take the column data to be located at the lower portion of non-negative axis, i.e. Rc=
Ri(l:End, c), it finds out all extreme points (both ends boundary point is not counting extreme point) of Rc and finds first extreme point and be located at Rc
Coordinate value Rcd, Rcd be approximately BIMFiThe period of (x, y) main component in a column direction is arranged if Rc does not have extreme point
Rcd=1；
3.2)BIMFiThe global frequencies fR approximate calculation of (x, y) is：
In turn, extension fRw=[fRw fR].
4) k1 is determined to separate high fdrequency component and striped carrier component：
Work as k1When ≠ 0, high fdrequency component, which has determined that, to be finished, and next step is directly entered；
Work as k1When=0, if being less than 3 elements in Ep set, it can not determine whether minimum, return to step 2 to Ii+1
(x, y) continues to decompose；
If being more than or equal to 3 elements in Ep set, if Ep (1)<Ep (2), first energy are just minimum, then k1=i-
2, alternatively, if Ep (1)>Ep (2) and Ep (2)<Ep (3), can quantitative change dissolve existing minimum, then k1Otherwise=i-1 is returned to
Step 2 is to Ii+1(x, y) continues to decompose.
If 5) k1≠ 0 and i>k1+ 1, high fdrequency component, which has determined that, to be finished, and starts to determine k2With separate striped carrier component and
Otherwise low frequency component returns to step 2；Determine k2The step of be：
5.1) fR obtained with step 3.2 extends fRk, obtains fRk=[fRk fR], element number n_fRk；
5.2) if n_fRk ＜ 2, return to step 2；If n_fRk >=2,：
5.2.1) calculate BIMFi-1(x, y) and BIMFiThe global frequencies ratio of (x, y) is rf=fRk (i-1)/fRk (i), and
Extend rfk=[rfk rf]；Calculate BIMFiThe global amplitude of (x, y) isWherein E can be obtained by step 3, and be extended
ARk=[aRk aR] and enable its element number be n_aRk；Take the value for corresponding to identical BIMF in fRk with aRk, i.e. fRka=fRk
(n_fRk-n_aRk+1:End), and by aRk is corresponding with the element in fRka it is divided by, obtains amplitude-frequency and be combined into af=than collection
aRk/fRka；
5.2.2) if there is the element more than or equal to 2 in rfk, 5.2.3 is entered step；
If being not greater than the element equal to 2 in rfk, work as AvDistPreiWhen ＜ round [max (c, l)/3], step is returned to
2；Work as AvDistPreiWhen >=round [max (c, l)/3], if the value in af changes in monotone decreasing, k2=(k1+ 1)+1, it is no
Then, k2=i, BSEMD decompose stopping and k1、k2Confirmation work terminate, k=i；
5.2.3) if the element of af set is not more than 2, step 2 is returned to；
Otherwise：
As af (i-1)≤af (i-2) and af (i-1) ＜ af (i), i.e. af (i-1) is minimum, then k2=i-1, and
BSEMD decomposes stopping and k1、k2Confirmation work terminate, k=i；
As af (i-1) >=af (i-2) and af (i-1) ＞ af (i), i.e. af (i-1) is maximum：If ENERGY E p (i-1)
It is also a maximum, if AvDistPre at this timei＜ round [max (c, l)/3], then return to step 2；If AvDistPrei≥
Round [max (c, l)/3], then k2=i-1, BSEMD decompose stopping and k1、k2Confirmation work terminate, k=i；If ENERGY E p
(i-1) be not a maximum, then k2=i-2, BSEMD decompose stopping and k1、k2Confirmation work terminate, k=i；
As af (i-2) ≈ af (i-1) ≈ af (i), step 2 is returned to；Otherwise, if AvDistPrei≥round[max(c,
L)/3], BSEMD decomposes stopping and k1、k2Confirmation work terminate, k=i, simultaneously：If the value in af changes in monotone decreasing,
k2=n2_rf+ (k1+ 1), wherein n2_rf is that first element of the value more than 2 corresponds to the ordinal number that rf is concentrated in rf；If in af
Value be in monotone increasing, k2=i.
6) striped carrier component is post-processed, steps are as follows：
6.1) B (x, y)=BIMF is enabledi(x, y), i are initiated with (k1+1), revolve transformation (Hilbert with Hilbert
Spiral transform, HST) calculate B (x, y) two-dimentional instantaneous amplitude distribution, obtain sA (x, y)；With Otsu Threshold segmentations
SA (y, x) is divided into ii blocks, gained segmentation result to be denoted as oA (y, x) by method；Wherein, HST asks two-dimentional instantaneous amplitude can be by following formula
It obtains：
Wherein, F () and F-1() indicates two-dimensional Fourier transform and two-dimentional inverse Fourier transform respectively,For at the rotation phase equation of Fourier transform domain (Spiral phase function),For
Spectrum domain coordinate system；
6.2) operation is opened to oA (y, x) progress morphology and obtains openA (x, y), carrying out morphology to openA (x, y) closes
Operation obtains closeA (x, y), further, carries out morphological dilations to closeA (x, y) and obtains dilaA (x, y)；
6.3) all in detection dilaA (y, x) to mark the value for being and their coordinate, find B according to these coordinates
Analog value and they are set to 0 in (x, y)；
6.4) sig'(x, y are enabled)=sig (x, y)+B (x, y)；Enable ii=ii+1, and i=i+1；
If 6.5) i=k3(k3 set Ep (k1+1), Ep (k1+2) ..., Ep (k2) the corresponding subscript of maximum value), after
Processing terminates, and otherwise, returns to step 6.1, repeats the process；
6.6) it will be added with without the striped component by post-processing by post-processing,
Final enhanced striped carrier component PM=sig " (x, y).
Step 2：Hilbert transform is carried out to striped carrier component and obtains wrapped phase and instantaneous frequency fins(x, y) simultaneously
Quality Map Q (x, y) is established using instantaneous frequency, wherein instantaneous frequency establishes Quality Map according to following formula：
Q (x, y)=| fins(x,y)-f0|
F in formula0For striped carrier frequency.Instantaneous frequency is calculated as：PM components are pressed using Matlab tool functions instfreq
Row seeks one-dimensional instantaneous frequency, result fins1(x, y), since instantaneous frequency can not be corresponded with pixel, so do as
Lower simplification is to be mapped：
fins=[fins,fins(:,1:2)]
Striped carrier frequency estimates f0Mode is：Equal is not asked by the instantaneous frequency of the one-row pixels point of object modulation to bar graph
Value.Wrapped phase and the Quality Map difference of acquisition are as shown in Figure 4 and Figure 5.
Step 3：Edge detection is carried out to high fdrequency component, edge detection mode is Canny operator edge detections, is used
Matlab tool functions edge can be obtained, and edge pixel values are 1 in note two values matrix MASK1；Low frequency component is utilized
MATLAB tool function otsu multi-threshold segmentations, using matlab tool function otsu, multi-threshold is generally 5, and gray value is minimum
A kind of pixel be shadow region, shadow region pixel value is 1 in note two values matrix MASK2.As a result respectively such as Fig. 6, Fig. 7 institute
Show.
And then, the two merges, using matrix or operation, i.e.,：
MASK3=MASK1 | | MASK2
Then morphology is embodied as expansive working：
Structural element BI is circle, and radius size generally chooses 5 pixel sizes, low quality region two values matrix MASK
Intermediate value forms for 1 pixel.
Step 4：Each connected component labeling, the specific steps are：
7.1) each 4 connection in the region that the quick recognition matrix MASK values of MATLAB tool functions bwlabel are 0 is used
Region, in each region grey scale pixel value be respectively 1,2 ..., m (m is high quality areal)；
7.2) each 4 connection in the region that the quick recognition matrix MASK values of MATLAB tool functions bwlabel are 1 is used
Region, in each region grey scale pixel value be respectively m+1, m+2 ..., m+n (n is low quality areal)；
7.3) two matrixes of the above results and as each connected component labeling result (are indicated, MARK values with matrix MARK
It is 1,2 ..., m indicates that each high quality region, MARK values are m+1 respectively, and m+2 ... m+n indicate each low mass region respectively
Domain), the results are shown in Figure 8 (m=2 in this example, n=2)；
Construct the low quality region of extension：Include boundary pixel point, boundary pixel again by low quality region described in MARK
Point is uniformly selected in adjacent high quality region, is illustrated below：It is low quality region if Fig. 9 is white area, as Figure 10 is white
Color region is the low quality region of extension.
Phase flood sprawling expansion, is embodied as in each high quality region：
8.1) it is starting point to choose any point in high quality, and queue is added in its coordinate；
8.2) belong to the phase of the pixel in high quality region in four neighborhood of expansion queue heads element, expansion mode presses such as following table
Show that formula is unfolded；
Current required expansion phase, ψwrpCurrent package phase,Phase is unfolded in former point.
8.3) queue is added in the phase coordinate being unfolded in step 8.2), 8.1) queue heads element dequeue continues to execute
Until queue is sky.
The low quality region Quality Map method expansion of each extension, process description are：Quality maximum point is made in chosen area
For starting point, the phase for the point being not yet unfolded in four adjoint points in the region is unfolded and the information of these points is stored in an array
In, by the first information exchange in the point of mass value in array and array, first point in array is taken out into update array,
It repeats above operation, until there is no element in array.
Step 5：Two adjacent region phases are successively merged, until all areas phase all merges completion, wherein adjacent
Region phase merges algorithm and is expressed as with following formula：.
Round () is the operation that rounds up in formula,It is unfolded as a result, ψ in high quality region for boundary element phasejFor
Boundary element is unfolded in the low quality region of extension as a result, n is boundary pixel sum.Final result is as shown in figure 11.
Embodiment is merely illustrative of the invention's technical idea, and cannot limit protection scope of the present invention with this, it is every according to
Technological thought proposed by the present invention, any change done on the basis of technical solution, each falls within the scope of the present invention.
Claims (9)
1. a kind of subregional single width fast phase method of deploying, which is characterized in that include the following steps：
(1) it uses BSEMD to carry out multi-resolution decomposition to bar graph I (x, y), obtains the exploded view BIMF of k width different scalesi(x,
y)：
Wherein, BIMFiFor intrinsic mode function, r (x, y) is surplus, and determines parameter k1、k2To divide high fdrequency componentThe striped carrier component of enhancingAnd low frequency component
(2) Hilbert transform is carried out to striped carrier component and obtains wrapped phase and instantaneous frequency fins(x, y), and utilize wink
When frequency fins(x, y) establishes Quality Map Q (x, y)；
(3) edge detection is carried out to high fdrequency component, shadow region identification is carried out to low frequency component, both fusions are as a result, again to melting
Result applied morphology after conjunction operates to form low quality region, remaining region is high quality region；
(4) label distinguishes each connected domain, using flood sprawling expansion phase, the low quality of extension in each high quality region
Phase is unfolded using Quality Map method in region；
(5) adjacent high and low quality region phase merges, until ultimately forming complete continuous phase.
2. subregional single width fast phase method of deploying according to claim 1, which is characterized in that in step (2), press
Quality Map Q (x, y) is obtained according to following formula：
Q (x, y)=| fins(x,y)-f0|
Wherein, f0For striped carrier frequency.
3. subregional single width fast phase method of deploying according to claim 1, which is characterized in that in step (3), adopt
Edge detection is carried out to high fdrequency component with Canny operators, edge pixel values in two values matrix MASK1, which are set to 1 as edge, to be examined
Survey result.
4. subregional single width fast phase method of deploying according to claim 3, which is characterized in that right in step (3)
Low frequency component carries out shadow region and knows method for distinguishing：Multi-threshold segmentation is carried out using otsu methods to low frequency component, as a result middle gray scale
It is shadow region to be worth minimum pixel forming region, and shadow region pixel value in two values matrix MASK2, which is set to 1, is used as shadow region
Domain recognition result.
5. subregional single width fast phase method of deploying according to claim 4, which is characterized in that in step (3), institute
It refers to that two values matrix MASK1 and MASK2 progress or operation are obtained MASK3 to state fusion.
6. subregional single width fast phase method of deploying according to claim 4, which is characterized in that in step (3), institute
It is as follows to state morphological operation：
Wherein, BI is circular structural element, and the low quality region is made of the pixel that two values matrix MASK intermediate values are 1.
7. subregional single width fast phase method of deploying according to claim 6, which is characterized in that in step (4), mark
The step of remembering each connected domain is as follows：
(401) each 4 connected region in the region that the quick recognition matrix MASK values of MATLAB tool functions bwlabel are 0 is used
Domain, in each region grey scale pixel value be respectively 1,2 ..., m, wherein m are high quality areal；
(402) each 4 connected region in the region that the quick recognition matrix MASK values of MATLAB tool functions bwlabel are 1 is used
Domain, in each region grey scale pixel value be respectively m+1, m+2 ..., m+n, wherein n is low quality areal；
(403) result of step (401) and (402) is done into matrix and operation, operating result be each connected component labeling as a result,
It is indicated with matrix MARK, MARK values are 1,2 ..., and m indicates that each high quality region, MARK values are m+1, m+2 ..., m+n respectively
Each low quality region is indicated respectively.
8. subregional single width fast phase method of deploying according to claim 7, which is characterized in that in step (4), institute
State the meaning in the low quality region of extension：Low quality region includes boundary pixel point in matrix MARK, and boundary pixel point is unified
It is selected in adjacent high quality region.
9. subregional single width fast phase method of deploying according to claim 1, which is characterized in that in step (5), phase
Adjacent two region phase merging process are indicated with following formula：
Wherein, KjFor phase amalgamation result, round () is the operation that rounds up,It is boundary element phase in high quality region
Interior expansion is as a result, ψjIt is unfolded as a result, n is boundary pixel sum in the low quality region of extension for boundary element.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810352746.7A CN108680119B (en) | 2018-04-19 | 2018-04-19 | Partitioned single-amplitude fast phase unwrapping method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810352746.7A CN108680119B (en) | 2018-04-19 | 2018-04-19 | Partitioned single-amplitude fast phase unwrapping method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108680119A true CN108680119A (en) | 2018-10-19 |
| CN108680119B CN108680119B (en) | 2020-06-12 |
Family
ID=63801468
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810352746.7A Active CN108680119B (en) | 2018-04-19 | 2018-04-19 | Partitioned single-amplitude fast phase unwrapping method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108680119B (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109472834A (en) * | 2018-10-23 | 2019-03-15 | 桂林电子科技大学 | A kind of Kalman filtering phase developing method based on wavelet transformation |
| CN112923847A (en) * | 2021-01-21 | 2021-06-08 | 广东工业大学 | Local sine auxiliary grating ruler measurement error adaptive compensation method |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2198780A2 (en) * | 2008-12-19 | 2010-06-23 | Sirona Dental Systems GmbH | Method and device for optical scanning of three-dimensional objects by means of a dental 3D camera using a triangulation method |
| CN102322822A (en) * | 2011-08-08 | 2012-01-18 | 西安交通大学 | Three-dimensional measurement method for triple-frequency color fringe projection |
| US20150062558A1 (en) * | 2013-09-05 | 2015-03-05 | Texas Instruments Incorporated | Time-of-Flight (TOF) Assisted Structured Light Imaging |
| CN105354877A (en) * | 2015-11-09 | 2016-02-24 | 北京航空航天大学 | Three-dimensional grid processing method based on empirical mode decomposition and Hilbert spectrum calculation of space filling curve |
| CN106170679A (en) * | 2015-09-02 | 2016-11-30 | 深圳大学 | A kind of phase error compensation method and device |
-
2018
- 2018-04-19 CN CN201810352746.7A patent/CN108680119B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2198780A2 (en) * | 2008-12-19 | 2010-06-23 | Sirona Dental Systems GmbH | Method and device for optical scanning of three-dimensional objects by means of a dental 3D camera using a triangulation method |
| CN102322822A (en) * | 2011-08-08 | 2012-01-18 | 西安交通大学 | Three-dimensional measurement method for triple-frequency color fringe projection |
| US20150062558A1 (en) * | 2013-09-05 | 2015-03-05 | Texas Instruments Incorporated | Time-of-Flight (TOF) Assisted Structured Light Imaging |
| CN106170679A (en) * | 2015-09-02 | 2016-11-30 | 深圳大学 | A kind of phase error compensation method and device |
| CN105354877A (en) * | 2015-11-09 | 2016-02-24 | 北京航空航天大学 | Three-dimensional grid processing method based on empirical mode decomposition and Hilbert spectrum calculation of space filling curve |
Non-Patent Citations (2)
| Title |
|---|
| FEIPENG DA ET.: "A fast, accurate phase unwrapping method for wavelet-transform profilometry", 《OPTICS COMMUNICATIONS》 * |
| 梁曼等: "瞬时频率引导的小波变换轮廓术相位解包裹技术", 《光电子·激光》 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109472834A (en) * | 2018-10-23 | 2019-03-15 | 桂林电子科技大学 | A kind of Kalman filtering phase developing method based on wavelet transformation |
| CN109472834B (en) * | 2018-10-23 | 2023-04-14 | 桂林电子科技大学 | Kalman filtering phase expansion method based on wavelet transformation |
| CN112923847A (en) * | 2021-01-21 | 2021-06-08 | 广东工业大学 | Local sine auxiliary grating ruler measurement error adaptive compensation method |
| CN112923847B (en) * | 2021-01-21 | 2022-03-18 | 广东工业大学 | Local sine auxiliary grating ruler measurement error adaptive compensation method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108680119B (en) | 2020-06-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US10325151B1 (en) | Method of extracting image of port wharf through multispectral interpretation | |
| CN101853333B (en) | Method for picking marks in medical robot navigation positioning images | |
| CN103729251B (en) | Concurrent computation optical bar chart phase extraction method | |
| CN103414861B (en) | A kind of method of projector frame self-adaptive Geometry rectification | |
| CN106846344A (en) | A kind of image segmentation optimal identification method based on the complete degree in edge | |
| CN102521859B (en) | Reality augmenting method and device on basis of artificial targets | |
| CN106709947A (en) | RGBD camera-based three-dimensional human body rapid modeling system | |
| CN109087345A (en) | Pallet recognition methods and automated guided vehicle based on ToF imaging system | |
| CN109146838A (en) | A kind of aobvious band adhering chromosome dividing method of the G merged based on geometrical characteristic with region | |
| CN108680119A (en) | A kind of subregional single width fast phase method of deploying | |
| CN110207620B (en) | Three-dimensional reconstruction method for determining optical series of digital grating projection structure through different frequencies | |
| CN105184251B (en) | A kind of recognition methods in the wawter bloom region based on high score No.1 satellite image and device | |
| US20070248268A1 (en) | Moment based method for feature indentification in digital images | |
| CN105547190B (en) | 3 D measuring method and device based on double angle unifrequency fringe projections | |
| CN102930279B (en) | For the image-recognizing method that product quantity detects | |
| CN110853081A (en) | Ground and airborne LiDAR point cloud registration method based on single-tree segmentation | |
| CN108664733A (en) | Seamed edge characteristic surface topology approaches method for reconstructing | |
| CN109522904B (en) | Rule farmland extraction method based on remote sensing data | |
| CN111881801B (en) | Newly-added construction land remote sensing monitoring method and equipment based on invariant detection strategy | |
| CN105654423A (en) | Area-based remote sensing image registration method | |
| CN111932632A (en) | Phase correction method in three-dimensional reconstruction of mechanical part | |
| CN105588518B (en) | Three-dimensional appearance acquisition methods based on double angle multi-frequency fringe projections and device | |
| CN107918938A (en) | A kind of matching process and device of point set and point set | |
| CN106931905A (en) | A kind of digital Moiré patterns phase extraction method based on nonlinear optimization | |
| Da et al. | A fast, accurate phase unwrapping method for wavelet-transform profilometry |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |