CN108680119B

CN108680119B - Partitioned single-amplitude fast phase unwrapping method

Info

Publication number: CN108680119B
Application number: CN201810352746.7A
Authority: CN
Inventors: 王辰星; 聂文祥; 达飞鹏
Original assignee: Southeast University
Current assignee: Southeast University
Priority date: 2018-04-19
Filing date: 2018-04-19
Publication date: 2020-06-12
Anticipated expiration: 2038-04-19
Also published as: CN108680119A

Abstract

The invention discloses a partitioned single-amplitude fast phase unwrapping method. Firstly, performing self-adaptive multi-scale decomposition on a fringe pattern by using a two-dimensional sine-assisted empirical mode decomposition method (BSEMD), wherein an obtained high-frequency result component is used for detecting the edge of an object in a fringe, an obtained intermediate fringe carrier frequency component obtains a wrapping phase and instantaneous frequency through Hilbert transformation, and an obtained low-frequency component is used for recognizing a shadow region of the fringe pattern; then fusing the edge and shadow regions and performing morphological operation to form a low-quality region; then, expanding a high-quality area by using a flood spreading mode, and establishing a quality diagram by using instantaneous frequency to guide the phase expansion of the expanded low-quality area; finally, the phases of all the areas are combined to form a complete continuous phase. The invention realizes the rapid unfolding of the wrapping phase by applying different phase unfolding strategies in the divided areas, has high efficiency and high precision in the whole method, and has larger application prospect in a single stripe three-dimensional rapid measurement system.

Description

Partitioned single-amplitude fast phase unwrapping method

Technical Field

The invention belongs to the technical field of three-dimensional information reconstruction, and particularly relates to a partitioned single-amplitude rapid phase unwrapping method.

Background

The fringe projection measurement technology has become a popular three-dimensional surface shape measurement technology due to the advantages of non-contact, rapidness, accuracy and the like, and is widely applied to the fields of national defense safety, advanced manufacturing, reverse engineering and cultural relic protection. The optical measurement technology of a plurality of stripes can realize high-precision measurement with less condition constraint. However, this method needs to keep the object still to project multiple images, and cannot realize dynamic three-dimensional measurement, so a single fringe projection technique is produced.

The fringe projection technique projects sinusoidal grating fringes toward a measured object and photographs a fringe pattern modulated by the object surface, three-dimensional information of the object surface is contained in the fringe phase, and thus demodulation of the fringe image is required. The phase acquisition technology based on single fringe projection mainly comprises the following steps: fourier transform, windowed fourier transform, S-transform, wavelet transform, and two-dimensional empirical mode decomposition (BSEMD), among others. The BSEMD is used as a new self-adaptive time-frequency analysis method, and has better time-frequency analysis capability than other methods due to being driven by self signals when analyzing unsteady-state signals, and the application scene of the invention is based on the BSEMD.

These transform domain based phase demodulation are all obtained by arctan transform, so that only the dominant phase is obtained in the range of-pi, pi), and therefore needs to be unwrapped to continuous phase. The phase unwrapping is generally accomplished by only comparing whether the wrapped phase difference between two adjacent points exceeds pi to compensate for the 2 pi offset. However, the conventional line-by-line deployment method is prone to the "pull-off" phenomenon due to poor phasing caused by undersampling, noise, invalid shadow areas, sharp object surfaces, and the like.

The quality map method is the most common one of the phase unwrapping methods, and the quality map describes the phase quality of each pixel point in the fringe pattern so as to guide the phase integral unwrapping. Although the path unfolding method guided by the quality map can control the propagation of phase errors, the quality of each point is compared to select high-quality points to be preferentially unfolded, so that the running time of an algorithm is too long, and the rapid characteristic of single fringe projection is greatly reduced.

Disclosure of Invention

In order to solve the technical problems of the background art, the invention aims to provide a partitioned single-frame fast phase unwrapping method, and the operating time of a quality map method is reduced by unfolding phases in a partitioned mode.

In order to achieve the technical purpose, the technical scheme of the invention is as follows:

a partitioned single-amplitude fast phase unwrapping method includes the following steps:

(1) carrying out multiscale decomposition on the fringe image I (x, y) by using BSEMD to obtain k exploded views BIMF with different scales_i(x,y)：

Wherein, BIMF_iFor the natural mode function, r (x, y) is the margin, and the parameter k is determined₁、k₂To divide the high frequency component

Enhanced fringe carrier frequency component

And low frequency components

(2) Hilbert transform is performed on the striated carrier frequency component to obtain a wrapped phase and an instantaneous frequency f_ins(x, y) and using the instantaneous frequency f_ins(x, y) establishing a quality map Q (x, y);

(3) performing edge detection on the high-frequency component, performing shadow region identification on the low-frequency component, fusing the two results, applying morphological operation to the fused result to form a low-quality region, and taking the rest regions as high-quality regions;

(4) marking and distinguishing each connected domain, adopting flood spreading and expanding phases in each high-quality area, and adopting a quality map method to expand phases in the expanded low-quality area;

(5) the phases of adjacent high and low quality regions are combined until a complete continuous phase is finally formed.

Further, in the step (2), a mass chart Q (x, y) is obtained according to the following formula:

Q(x,y)＝|f_ins(x,y)-f₀|

wherein f is₀Is the fringe carrier frequency.

Further, in step (3), edge detection is performed on the high-frequency component by using a Canny operator, and an edge pixel value in the binary matrix MASK1 is set to 1 as an edge detection result.

Further, in step (3), the method for identifying the shadow region of the low-frequency component comprises the following steps: and performing multi-threshold segmentation on the low-frequency components by using an otsu method, wherein a pixel forming area with the lowest gray value in the result is a shadow area, and setting the pixel value of the shadow area in the binary matrix MASK2 to be 1 as a shadow area identification result.

Further, in the step (3), the fusion refers to performing or operation on the binary matrix MASK1 and MASK2 to obtain MASK 3.

Further, in step (3), the morphological operations are as follows:

wherein, BI is a circular structural element, and the low-quality area is composed of pixels with a median value of 1 in a binary matrix MASK.

Further, in step (4), the step of labeling each connected domain is as follows:

(401) rapidly identifying each 4-connected region in the region with the matrix MASK value of 0 by using an MATLAB tool function bwleael, wherein the gray value of pixels in each region is 1,2, … and m, and m is the number of high-quality regions;

(402) rapidly identifying each 4-connected region in the region with the matrix MASK value of 1 by using an MATLAB tool function bwleabel, wherein the gray values of pixels in each region are m +1, m +2, … and m + n respectively, and n is the number of the low-quality regions;

(403) and (3) performing matrix and operation on the results of the steps (401) and (402), wherein the operation result is the labeling result of each connected domain, the labeling result is represented by a matrix MARK, the MARK values are 1,2 and …, m respectively represents each high-quality region, and the MARK values are m +1, m +2, …, and m + n respectively represents each low-quality region.

Further, in step (4), the meaning of the extended low-quality region: the low-quality area in the matrix MARK comprises boundary pixel points, and the boundary pixel points are uniformly selected in the adjacent high-quality areas.

Further, in step (5), the phase merging process of two adjacent regions is represented by the following formula:

wherein, K_jAs a result of the phase combining, round () is a rounding operation,

for the result of the spreading of the phase of the boundary element in the high-quality region,. psi_jN is the total number of boundary pixels as a result of the boundary elements being spread out within the expanded low-quality region.

Adopt the beneficial effect that above-mentioned technical scheme brought:

(1) the invention utilizes BSEMD to acquire the wrapping phase and utilizes the intermediate result instantaneous frequency to establish the quality diagram, thereby greatly reducing the time for establishing the quality diagram.

(2) The method adopts self-adaptive region generation, the region division depends on the surface characteristics and the illumination angle of the object, and compared with a fixed region division method, the method has more flexibility and more scientific and reasonable processing on the region boundary.

(3) Compared with global quality map method expansion, the method has higher expansion speed, so that the single stripe is more suitable for fast three-dimensional measurement.

Drawings

FIG. 1 is an overall flow chart of the present invention;

FIG. 2 is a stripe diagram of an embodiment;

FIG. 3 is a flow chart of the BSEMD algorithm;

FIG. 4 is a graph of example wrapped phase;

FIG. 5 is a quality chart established by the embodiment;

FIG. 6 is a schematic diagram of edge detection according to an embodiment;

FIG. 7 is a diagram illustrating the result of detecting a shadow area according to an embodiment;

FIG. 8 is a schematic illustration of an embodiment of area ID numbering;

FIGS. 9 and 10 are schematic diagrams for auxiliary illustration of boundary elements according to an embodiment;

FIG. 11 is a diagram showing the final phase unwrapping result of the embodiment.

Detailed Description

The technical scheme of the invention is explained in detail in the following with the accompanying drawings.

As shown in fig. 1, in the single-amplitude fast phase unwrapping method for segmenting a region, firstly, a two-dimensional sine-assisted empirical mode decomposition (BSEMD) is used to perform adaptive multi-scale unwrapping on a fringe pattern, an obtained high-frequency result component is used to detect the edge of an object in a fringe, an obtained intermediate fringe carrier frequency component obtains a wrapping phase and an instantaneous frequency through hilbert transform, and an obtained low-frequency component is used to identify a shadow region of the fringe pattern. The edges and shadow regions are then fused and morphologically manipulated to form low quality regions. And then, expanding the high-quality area by using a flood spreading mode, and establishing a quality map by using instantaneous frequency to guide the phase expansion of the low-quality area. Finally, the phases of all the areas are combined to form a complete continuous phase. The detailed procedure of the method is described below.

The method comprises the following steps: the high frequency component is divided from the stripe by decomposing the stripe map Im (x, y) (as shown in FIG. 2) by the BSEMD algorithm

Enhanced fringe carrier frequency component PM, low frequency component

The algorithm flow chart is shown in fig. 3, and the operation steps are described as follows:

1) initializing attribute quantity:

I_i(x, y) ═ Im (x, y), each exploded view BIMF_iThe ordinal number i of (x, y) is 1, c, l is the stripe width and height, and k1 is 0;

BIMF_istatistical average period AvDistPre ═ 1 for (x, y), each BIMF_i(x, y) the set Exe of statistical averaging periods is an empty set;

BIMF_ienergy E of (x, y) is 0, each BIMF_i(x, y) the energy set Ep is an empty set;

BIMF_iglobal frequency fR of (x, y) is 0, each BIMF_i(x, y) the global frequency set fRw is an empty set;

fRk is empty for storage

Then each BIMF_iA global frequency of (x, y);

aRk is empty for storage

Then each BIMF_iA global magnitude of (x, y);

rfk is empty for storage

Then the global frequency ratio of two adjacent BIMF;

the final fringe carrier component is sig-0.

2) Using BSEMDFair I_i(x, y) performing a single-layer scale decomposition, wherein the BSEMD single-layer decomposition function has an input of I_i(x, y) and a statistical average period parameter AvDistPree, the output is BIMF_i(x, y), the BSEMD algorithm specifically comprises the following steps:

2.1) to I_i(x, y) performing single-layer decomposition by using an Enhanced Fast Empirical Mode Decomposition (EFEMD) method, wherein the input of an EFEMD single-layer decomposition function is I_i(x, y) and a statistical average period parameter AvDistPree, the output is BIMF^out(x, y) and a new statistical averaging period, avdistpor 0; is provided with

Is I_i(x, y) initial exploded view, in particular, the single layer decomposition steps of EFEMD are as follows:

2.1.1) pairs I_i(x, y) mirror extension to produce I_e(x, y); specifically, the Matlab-based mirror extension steps are as follows:

determining the continuation size as es ═ round [ min (c, l)/10 ];

extracting an outer frame with the width of es of the image Im (x, y), and solving a mirror surface of the outer frame as a peripheral extension value, namely: horizontal top extension value ext1 ═ flipud [ Im (1: es,1: c)]Wherein the flit (-) is a matrix up-down flip function carried in Matlab; horizontal bottom extension value ext2 ═ flash [ Im (l-es +1: l,1: c)](ii) a The left continuation value is ext3 ═ fliplr [ Im (1: l,1: es)]Wherein the fliplr (-) is a matrix left-right turning function carried in Matlab; the right extension is ext4 ═ fliplr [ Im (1: l, c-es +1: c)](ii) a The upper left continuation value is ext5 ═ fliplr { rot90[ Im (1: es )]R, where rot90 (-) is a Matlab self-carrying function to realize that the matrix rotates 90 degrees counterclockwise; the upper right continuation value is ext6 ═ Im (1: es, c-es +1: c)^T，(·)^TTransposing the matrix; lower left corner extension value ext7 ═ Im (l-es +1: l,1: es)^T(ii) a The lower right continuation value is ext8 ═ fliplr { rot90[ Im (l-es +1: l, c-es +1: c)]}；

Combining all the extension values of the above steps with Im (x, y), the final result after mirror extension is:

2.1.2) pairs I_e(x, y) two morphological dilation operations were performed separately, namely:

MaxD(x,y)＝imdilate(I_e(x,y),ms)

MinD(x,y)＝-imdilate(-I_e(x,y),ms)

where ms is strel ('disk', AvDistPre), strel (·) is a self-contained function in the Matlab toolbox to create a structural element of a specified shape, and imdilate (·) is a self-contained morphological expansion function in Matlab; to obtain I_eThe maximum distribution of (x, y) is:

MaxM(x,y)＝～(I_e(x,y)-MaxD(x,y))

I_ethe distribution of minimum values of (y, x) is:

MinM(x,y)＝～(I_e(x,y)-MinD(x,y))

wherein, the 'to' is 'logical not' operation, namely, the value of 0 in the matrix is set as 1, and the value of 0 is set as 0;

2.1.3) calculate the average spatial distance of all extrema:

wherein round (·) is an integer fetch operation; if the AvDistNew is larger than the AvDistNew, the AvDistNew is made to be AvDistNew, otherwise, the AvDistNew is made to be 2 multiplied by the AvDistNew;

2.1.4) averaging the two morphological dilation results of step 2.1.2 and performing a two-dimensional convolution smoothing operation on the obtained average:

where mC ═ fspecial ("disk", round (AvDistPre/2)), conv2(·) is a two-dimensional convolution operation function in the Matlab toolbox, and fspecial (·) is a function that establishes a predefined filter operator in the Matlab toolbox.

2.1.5) step (b)Smoothing results obtained in step 2.1.4 from I_e(x, y) minus:

B(x,y)＝I_e(x,y)-SM(x,y)

the values of the continuation region in B (x, y) are removed according to step 2.1.1, resulting in a BIMF^out(x, y) is an exploded view of a layer obtained using an EFEMD which ultimately outputs a BIMF^out(x, y) and AvDistPre0, note

Is an initial exploded view;

2.2) obtained according to step 2.1

The two-dimensional sinusoidal distribution s (x, y) is designed as:

s(x,y)＝acos(2πfy)·cos(2πfx)

wherein the amplitude a is designed as the maximum value of all elements of the two-dimensional matrix of the exploded view:

max (-) is the max operation; the frequency f can be obtained by the following formula:

2.3) reacting s (x, y) with each other

Add and subtract, i.e.:

to pair

Performing EFEMD single-layer decomposition with input parameter of AvDistPree to obtain a first intermediate exploded view

Simultaneously obtaining the output new parameter AvDistPree^s1(ii) a Likewise, pair

EFEMD single layer decomposition (input parameter is still AvDistPree) is performed to obtain a second intermediate exploded view

And obtain newly output AvDistPree^s2；

2.4) to the resultant

And

and (3) averaging:

BIMF_i(x, y) is the final exploded view of the image Im (x, y); then, for two output parameters AvDistPree^s1And AvDistPree^s2Averaging, and making:

wherein round (·) is integer operation, and further, expand Exe is Exe ═ Exe AvDistPree_i]For AvDistPree_iThe following modifications were made:

(if AvDistPree)_i＜2Exe(end))

Wherein Exe (end) refers to the last value in the Exe set; mixing BIMF_i(x, y) from I_i(x, y) subtracting, let: i is_i+1(x,y)＝I_i(x,y)-BIMF_i(x,y)，I_i+1(x, y) as input graph for next decompositionLike this.

3) Calculating the BIMF_i(x, y) two-dimensional autocorrelation function:

wherein (tau)₁,τ₂) The time delay of (x, y) corresponding to the direction when the autocorrelation function is obtained is more than or equal to 0 and is more than or equal to tau₁≤2x-1，0≤τ₂≤2y-1；

According to R_i(τ₁,τ₂) Estimate the BIMF_iThe energy value of (x, y) is approximately:

E＝max[R_i(τ₁,τ₂)]

max [. cndot ] is the max operation; ep is extended to Ep ═ Ep E ];

according to R_i(τ₁,τ₂) Estimate the BIMF_iGlobal frequency fR of (x, y) as follows:

3.1) extraction of R_i(τ₁,τ₂) And the row data is taken to be located in the right half of the non-negative axis, i.e., Rl-R_i(l, c: end), calculating all extreme points of Rl (boundary points at two ends do not calculate the extreme points) and finding out a coordinate value Rld of which the first extreme point is positioned at Rl, wherein Rld is approximate to BIMF_i(x, y) a period of a main component in a row direction, and if Rl has no extreme point, setting Rld to 1;

similarly, R is extracted_i(τ₁,τ₂) And taking the lower half of the row of data on the non-negative axis, i.e. Rc ═ R_i(l: end, c), calculating all extreme points of Rc (boundary points at two ends do not calculate the extreme points) and finding out coordinate value Rcd of the first extreme point at Rc, wherein Rcd is approximate to BIMF_i(x, y) the period of the main component in the column direction, and if Rc has no extreme point, setting Rcd to 1;

3.2)BIMF_ithe global frequency fR of (x, y) is approximately calculated as:

further, extension fRw ═ fRw fR.

4) K1 is determined to separate the high frequency component from the fringe carrier frequency component:

when k is₁When the high-frequency component is determined to be finished, directly entering the next step;

when k is₁If the number of the elements in the Ep set is less than 3 when the number of the elements is 0, the minimum value cannot be determined, and the step 2 returns to the step I_i+1(x, y) continuing the decomposition;

if Ep is greater than or equal to 3 elements in the set, if Ep (1)<Ep (2), the first energy is extremely small, k₁I-2, or, if Ep (1)>Ep (2) and Ep (2)<Ep (3), i.e., the energy change, occurs to a minimum value, k₁Otherwise, go back to step 2 for I_i+1(x, y) the decomposition is continued.

5) If k is₁Not equal to 0 and i>k₁+1, high frequency component has already been confirmed, begin to confirm k₂Separating the frequency component and the low-frequency component by the stripe, otherwise, returning to the step 2; determining k₂Comprises the following steps:

5.1) extension fRk with fR obtained in step 3.2 to give fRk ═ fRk fR, with the number of elements n _ fRk;

5.2) if n _ fRk is less than 2, returning to the step 2; if n _ fRk ≧ 2, then:

5.2.1) calculating the BIMF_i-1(x, y) and BIMF_iThe global frequency ratio of (x, y) is rf fRk (i-1)/frk (i) and extends rfk ═ rfk rf](ii) a Calculating the BIMF_iA global amplitude of (x, y) of

Where E is available from step 3 and extended aRk ═ aR]And the number of elements is n _ aRk; taking fRk the value corresponding to aRk with the same BIMF, i.e. fRka is fRk (n _ fRk-n _ aRk +1: end), and dividing aRk with the element correspondence in fRka to obtain the amplitude-frequency ratio set af is aRk/fRka;

5.2.2) if element greater than or equal to 2 is present in rfk, go to step 5.2.3;

if rfk has no element greater than or equal to 2, then the AvDistPree_i＜round[max(c,l)/3]If yes, returning to the step 2; when AvDistPree_i≥round[max(c,l)/3]When k is monotonically decreasing if the value in af changes₂＝(k₁+1) +1, otherwise, k₂I, BSEMD decomposition stops and k₁、k₂The confirmation operation of (1) is finished, k is i;

5.2.3) if the element of the af set is not more than 2, returning to the step 2;

otherwise:

when af (i-1) ≦ af (i-2) and af (i-1) < af (i), i.e., af (i-1) is minimal, then k is₂I-1, and BSEMD decomposition stops and k₁、k₂The confirmation operation of (1) is finished, k is i;

when af (i-1) ≧ af (i-2) and af (i-1) > af (i), i.e., af (i-1) is maximum: if the energy Ep (i-1) is also a maximum, then if AvDistPree_i＜round[max(c,l)/3]Returning to the step 2; if AvDistPree_i≥round[max(c,l)/3]Then k is₂I-1, BSEMD decomposition stops and k₁、k₂The confirmation operation of (1) is finished, k is i; if the energy Ep (i-1) is not a maximum, k₂I-2, BSEMD decomposition stops and k₁、k₂The confirmation operation of (1) is finished, k is i;

when af (i-2) ≈ af (i-1) ≈ af (i), returning to step 2; otherwise, if AvDistPree_i≥round[max(c,l)/3]BSEMD decomposition stopped and k₁、k₂The confirmation operation of (2) is finished, and k is equal to i, and at the same time: if the value in af changes monotonically, then k₂＝n2_rf+(k₁+1), wherein n2_ rf is the ordinal number in the rf set corresponding to the first value greater than 2 in rf; if the value in af is monotonically increasing, k₂＝i。

6) Post-processing the carrier frequency component of the stripes, comprising the following steps:

6.1) let B (x, y) ═ BIMF_i(x, y), i starting with (k1+1), calculating the two-dimensional instantaneous amplitude distribution of B (x, y) using Hilbert Spiral Transform (HST) to obtain sA (x, y); dividing sA (y, x) into ii blocks by Otsu threshold segmentation method, and recording the division result as oA (y, x); wherein, the two-dimensional instantaneous amplitude value obtained by the HST can be obtained by the following formula:

wherein, F (-) and F^-1(. cndot.) denotes a two-dimensional Fourier transform and a two-dimensional inverse Fourier transform,

to be a Spiral phase function in the fourier transform domain,

is a frequency spectrum domain coordinate system;

6.2) carrying out morphological opening operation on the oA (y, x) to obtain openA (x, y), carrying out morphological closing operation on the openA (x, y) to obtain closeA (x, y), and further carrying out morphological expansion on the closeA (x, y) to obtain dilaA (x, y);

6.3) detecting all the values marked as 1 in the dilaA (y, x) and the coordinates thereof, finding the corresponding values in the B (x, y) according to the coordinates and setting the values as 0;

6.4) let sig' (x, y) be sig (x, y) + B (x, y); let ii be ii +1 and i be i + 1;

6.5) if i ═ k₃(k3 sets { Ep (k1+1), Ep (k1+2), …, Ep (k2) } subscript corresponding to the maximum value), the post-processing is finished, otherwise, the process returns to step 6.1, and the process is repeated;

6.6) adding the strip components which are subjected to post-processing and are not subjected to post-processing,

the final enhanced fringe carrier frequency component PM sig "(x, y).

Step two: hilbert transform is performed on the striated carrier frequency component to obtain a wrapped phase and an instantaneous frequency f_ins(x, y) and creating a quality map Q (x, y) using the instantaneous frequency, wherein the instantaneous frequency creates the quality map according to the following equation:

Q(x,y)＝|f_ins(x,y)-f₀|

in the formula f₀Is the fringe carrier frequency. The instantaneous frequency is calculated as: one-dimensional instantaneous frequency is obtained by row for the PM component by using a Matlab tool function instfreq, and the result is f_ins1(x, y), since the instantaneous frequency cannot correspond to the pixel point one-to-one, the following simplification is made to correspond:

f_ins＝[f_ins,f_ins(:,1:2)]

fringe carrier frequency estimation f₀The method comprises the following steps: and averaging the instantaneous frequency of a row of pixel points of the fringe pattern which are not modulated by the object. The acquired wrapped phase and mass maps are shown in fig. 4 and 5, respectively.

Step three: performing edge detection on the high-frequency component, wherein the edge detection mode is Canny operator edge detection and can be obtained by using a Matlab tool function edge, and the edge pixel value in a binary matrix MASK1 is recorded as 1; and (3) carrying out multi-threshold segmentation on the low-frequency components by using an MATLAB tool function otsu, wherein the multi-threshold is generally 5, the class of pixels with the lowest gray value is a shadow region, and the pixel value of the shadow region in a binary matrix MASK2 is recorded as 1. The results are shown in FIGS. 6 and 7, respectively.

Next, the two are fused, using a matrix or operation, i.e.:

MASK3＝MASK1||MASK2

the morphology was then embodied as a swelling operation:

the structural element BI is circular, the radius is generally 5 pixels, and the low-quality area is formed by using a binary matrix MASK with the median value of 1 pixel.

Step four: each communicated area is marked, and the method comprises the following specific steps:

7.1) rapidly identifying each 4-connected region in the region with the matrix MASK value of 0 by using a MATLAB tool function bwlebel, wherein the gray value of the pixel in each region is 1,2, … and m (m is the number of high-quality regions);

7.2) rapidly identifying each 4-connected region in the region with the matrix MASK value of 1 by using a MATLAB tool function bwlebel, wherein the gray values of pixels in each region are respectively m +1, m +2, …, and m + n (n is the number of low-quality regions);

7.3) the two matrix sums of the above results are the respective connected domain labeling results (denoted by matrix MARK with MARK values of 1,2, …, m respectively representing the respective high quality regions and m +1, m +2, … m + n respectively representing the respective low quality regions), and the results are shown in fig. 8 (in this example, m is 2 and n is 2);

constructing an extended low-quality region: the low-quality area in the MARK further includes boundary pixel points, and the boundary pixel points are uniformly selected in adjacent high-quality areas, as illustrated below: the white area is a low quality area as shown in fig. 9, and the white area is an extended low quality area as shown in fig. 10.

The phase flood spreading in each high-quality area is implemented as follows:

8.1) selecting any point in high quality as a starting point, and adding the coordinates of the point into a queue;

8.2) expanding the phases of the pixels belonging to the high-quality area in the four neighborhoods of the queue head element, wherein the expansion mode is expanded according to the following expression;

the currently sought unwrapped phase, psi_wrpThe current wrapped-phase is the phase of the current wrap,

the previous point unwrapps the phase.

8.3) adding the phase coordinates expanded in the step 8.2) into the queue, dequeuing the queue head element, and continuing to execute the step 8.1) until the queue is empty.

The quality map method of each expanded low-quality region is expanded, and the process is described as follows: selecting a point with the maximum quality in the area as a starting point, expanding the phase of the points which are not expanded in the four adjacent points in the area, storing the information of the points into an array, exchanging the information of the point with the first point in the array of the quality value in the array, taking out the first point in the array to update the array, and repeating the operation until no element exists in the array.

Step five: sequentially combining two adjacent region phases until all the region phases are combined, wherein the adjacent region phase combination algorithm is expressed by the following formula: .

Where round () is a rounding operation,

for the result of the spreading of the phase of the boundary element in the high-quality region,. psi_jN is the total number of boundary pixels as a result of the boundary elements being spread out within the expanded low-quality region. The final result is shown in fig. 11.

The embodiments are only for illustrating the technical idea of the present invention, and the technical idea of the present invention is not limited thereto, and any modifications made on the basis of the technical scheme according to the technical idea of the present invention fall within the scope of the present invention.

Claims

1. A partitioned single-amplitude fast phase unwrapping method is characterized by comprising the following steps:

Enhanced fringe carrier frequency component

And low frequency components

2. The partitioned single fast phase unwrapping method according to claim 1, wherein in step (2), the quality map Q (x, y) is obtained as follows:

Q(x,y)＝|f_ins(x,y)-f₀|

wherein f is₀Is the fringe carrier frequency.

3. The partitioned single fast phase unwrapping method according to claim 1, wherein in the step (3), Canny operator is adopted to perform edge detection on the high frequency component, and the edge pixel value in the binary matrix MASK1 is set to 1 as the edge detection result.

4. The partitioned single fast phase unwrapping method according to claim 3, wherein in the step (3), the shadow region identification is performed on the low frequency component by: and performing multi-threshold segmentation on the low-frequency components by using an otsu method, wherein a pixel forming area with the lowest gray value in the result is a shadow area, and setting the pixel value of the shadow area in the binary matrix MASK2 to be 1 as a shadow area identification result.

5. The partitioned single-frame fast phase unwrapping method according to claim 4, wherein in the step (3), the fusing is performed by OR-ing the binary matrix MASK1 with MASK2 to obtain MASK 3.

6. The partitioned single-frame fast phase unwrapping method according to claim 5, wherein in step (3), said morphological operations are as follows:

7. The partitioned single-frame fast phase unwrapping method according to claim 6, wherein in the step (4), the step of labeling each connected component is as follows:

8. The partitioned single fast phase unwrapping method according to claim 7, wherein in step (4), the meaning of said expanded low-quality region is: the low-quality area in the matrix MARK comprises boundary pixel points, and the boundary pixel points are uniformly selected in the adjacent high-quality areas.

9. The partitioned single fast phase unwrapping method according to claim 1, wherein in step (5), the phase merging procedure of two adjacent partitions is represented by the following formula: