当前位置:首页 >> 信息与通信 >>

JM Huntley


Unwrapping noisy phase maps by use of a minimum-cost-matching algorithm
J. R. Buckland, J. M. Huntley, and S. R. E. Turner

An algorithm for unwrapping noisy phase maps by mean

s of branch cuts has been proposed recently. These cuts join discontinuity sources that mark the beginning or end of a 2p phase discontinuity. After the placement of branch cuts, the unwrapped phase map is unique and independent of the unwrapping route. We show how a minimum-cost-matching graph-theory method can be used to ?nd the set of cuts that has the global minimum of total cut length, in time approximately proportional to the square of the number of sources. The method enables one to unwrap un?ltered speckle-interferometry phase maps at higher source densities 10.1 sources pixel212 than any previous branch-cut placement algorithm.

1.

Introduction

Analysis of fringes produced by optical interferometry is used to measure a range of physical parameters, for example, strain, displacement, and refractive index. Techniques such as moire ? , speckle, and holographic interferometry give a two-dimensional fringe pattern. The parameter of interest is related to the underlying phase distribution that is encoded in the fringe pattern. There are two main methods of extracting this distribution: 1i2 shifting the fringes through known phase increments and 1ii2 Fourier transformation of a single pattern containing carrier fringes.1 However, both methods result in a phase map wrapped onto the range 2p to p. This phase map must then be unwrapped 1an unknown multiple of 2p added to each pixel2 before it can be used to provide quantitative information. Figure 1 shows a wrapped phase map representing phases of 2p and 1p with black and white, respectively. The underlying 1unwrapped2 phase map is a vertically increasing uniform ramp. Noise has resulted in errors in the center of the phase map, causing a break in one of the 2p phase discontinuities at points 1 and 2.

Reference 2 gives a useful overview of phaseunwrapping algorithms developed so far. The simplest method involves comparing the phase of a pixel with that of one of its neighbors. If the fringe pattern is sampled at a sufficiently high spatial frequency 1at least twice the highest spatial frequency present in the pattern, in accordance with the Shannon sampling theorem2, then the true phase difference between neighboring pixels will be in the range 2p to p. A value outside this range indicates that a 2p discontinuity is crossed. Unwrapping from pixel i 2 1 to pixel i, we obtain the number of 2p phase discontinuities lying between the two pixels:

d1i2 5

3

F1i2 2 F1i 2 12 2p

4,

112

J. R. Buckland and J. M. Huntley are with the Cavendish Laboratory, University of Cambridge, Madingley Road, Cambridge CB3 OHE, UK; S. R. E. Turner is with the Department of Pure Mathematics and Mathematical Statistics, University of Cambridge, 16 Mill Lane, Cambridge CB2 1SB, UK. Received 17 August 1994; revised manuscript received 14 February 1995. 0003-6935@95@235100-09$06.00@0. r 1995 Optical Society of America. 5100 APPLIED OPTICS @ Vol. 34, No. 23 @ 10 August 1995

where square brackets denote rounding to the nearest integer and F1i2 is the phase of pixel i. d normally takes the value 0 3when F1i2 2 F1i 2 12 is in the range 2p to 1p4 but will take the value 11 or 21 when F1i2 2 F1i 2 12 is in the range p to 2p or 2p to 22p, respectively. If the phase is known at one pixel 1say, P in Fig. 12, one can calculate the phase at any other pixel in the image 1say, Q2 by following a path between the two pixels and counting the total number of discontinuities along the path:
N

n5

o d1i2.
i51

122

2pn is then subtracted from the phase of the Nth pixel. In the example shown in Fig. 1, na 5 23 for path a, but nb 5 22 for path b. The break in the

Fig. 1. Phase-unwrapping method. The phase at Q relative to P is obtained from the number of 2p discontinuities crossed by any path; noise results in discontinuity sources at points 1 and 2, causing unwrapping errors that are dependent on the path taken.

discontinuity caused by to the noise gives rise to an undesirable path dependence for the unwrapped phase and can cause the value at Q, which may be far from the noise, to be in error by a multiple of 2p. The branch-cut phase-unwrapping method3–5 involves placing a cut, which acts as a barrier to unwrapping, across the break in the 2p discontinuity, i.e., between points 1 and 2 in Fig. 1. All permitted unwrapping routes will then cross the same number of discontinuities and will result in the same unwrapped value. Points 1 and 2 can be regarded as discontinuity sources because they mark the beginning or end of lines of 2p discontinuity. They are easily detected by unwrapping counterclockwise around a closed loop of four pixels: n 5 61 if the loop encloses a source and zero otherwise. Points 1 and 2 are 11 and 21 sources, respectively. A fringe break that is due to noise will always result in such a 11@21 dipole pair. The formation of discontinuity sources can be seen in the example of Fig. 2. Figure 21a2 is a computergenerated cosine fringe pattern that contains two disks representing regions of low signal-to-noise ratio 1SNR2. Figure 21b2 shows the corresponding wrapped phase map calculated by the Fourier-transform method, and Fig. 21c2 shows the positions of the sources. The small disk results in the formation of a single dipole, whereas the larger one causes a cluster of ?ve dipoles. Each noise domain must have an equal number of positive and negative sources. One can see this by considering the two paths from P to Q in Fig. 21b2. These pass through noise-free regions of the image and so have the same value of n. The closed loop consists of one path from P to Q followed by the other path for P to Q; the two paths have opposite values of n, giving n 5 0 for the loop. It is easy to show 1see, for example, Ref. 42 that n for a closed loop equals the sum of the discontinuity sources enclosed by the loop, so that for every 11 source within the loop there must also be a 21 source. The correct distribution of cuts will therefore be a set of short lines linking 11 and 21 sources within each dipole pair. When a single noise domain contains

Fig. 2. 1a2 Computer-generated cosine fringe pattern with superimposed disks D1 and D2 to simulate regions of low SNR. 1b2 Wrapped phase map from the application of the Fourier-transform method to 1a2. 1c2 Source distribution for 1b2; the hollow and the ?lled circles are used to represent 11 and 21 sources, respectively.

more than one dipole 3e.g., that caused by D2 in Fig. 21b24, there may be several ways to join up the sources within the domain, all of which have approximately the same total cut length. However, although the choice of one particular set of cuts will affect the local phase distribution around the domain, the global phase map will be the same for all choices. Isolated monopoles can occur near the edge of the object if the partner source is off the ?eld of view, in which case the correct cut should be a short line from the monopole to the boundary. This approach will also cope with fringe patterns for which the underlying phase distribution is inherently discontinuous 1for example, objects containing cracks or holes, or ?uids with step changes in refractive index caused by shock waves2, provided boundary lines are set up to indicate these lines of true phase discontinuity. In general, the positions of the boundaries cannot be determined from the wrapped phase map alone, and additional information 1such as the modulation map or prior knowledge of the specimen geometry2 is required.
10 August 1995 @ Vol. 34, No. 23 @ APPLIED OPTICS 5101

The set up of the boundaries is problem speci?c and will not be addressed here; our aim in this paper is to deal reliably with the path-dependence problem once the boundaries are set up correctly. Reference 6 gives an overview of the methods that can be used to decide how to place the cuts, given an array of sources. All are based on the requirement that the branch cuts should be short; however, local search methods can result in a few long, physically unreasonable cuts. Although the phase map is now path independent, it contains 2p phase errors extending over the full ?eld of the image. In this paper we investigate the use of global minimization of the sum of the cut lengths, or the sum of the squares of the cut lengths, over the entire phase map as the criterion for placing the cuts. A theoretical justi?cation for minimizing the sum of squares of cut lengths in the case of speckle-interferometry phase maps is given in Ref. 7, based on experimental observations of source migration caused by speckle decorrelation. These criteria also seem intuitively appealing for phase maps obtained by other experimental techniques.
2. A. Description of the Algorithm Hungarian Algorithm

the optimal matching in the way described below, in which case the algorithm terminates. 152 If the maximum number of independent zeros in the matrix is L , N, then all the zeros in the matrix can be crossed out by L lines 1crossing out is discussed below2. After discarding any crossings from previous iterations, do this crossing out. 162 Find the minimum uncrossed entry, and subtract it from each uncrossed entry. 1This creates at least one more zero in the matrix2. Add the minimum uncrossed entry to each matrix element that is crossed in both directions. 172 Return to 142. A set of independent zeros in a matrix is a set of zero elements chosen so that no row or column contains more than one independent zero. A graph-theory method for ?nding a maximal set of independent zeros is discussed in Subsection 2.B. A crossing out is simply a horizontal or vertical line that crosses out the elements in a row or column. It can be shown that the above algorithm will always terminate. If there are L lines and D doublecrossed entries, then there are N 2 2 LN 1 D uncrossed entries. If L , N, this quantity is greater than D. Step 162 therefore involves more subtractions than additions so that the total of all the matrix entries decreases. Because all the entries are nonnegative integers, the algorithm must eventually terminate. That the algorithm does indeed provide an optimal matching can be proved by standard Lagrangian methods 1see, for example, Ref. 92. Once a maximal set of independent zeros has been found, the following algorithm is used to cross out all the zeros in the matrix 3step 152 above4: 112 For each row not containing an independent zero, cross out each zero vertically. 122 For each crossed-out independent zero, ?nd all the other uncrossed zeros in the same row. 132 If at least one uncrossed zero has been found in step 122, then cross it out vertically and go to 122. 142 Cross out all remaining uncrossed independent zeros horizontally. A matrix containing a maximal set of L independent zeros can have all the zeros crossed out by L lines, each passing through one independent zero. This fact, and the method of crossing out the zeros described above, can be viewed as special cases of the Ford–Fulkerson algorithm.9 Given that such a crossing out exists, a brief explanation of steps 112–142 is as follows. A zero in a row that contains no independent zeros cannot be crossed out horizontally because the line would not pass through an independent zero. A nonindependent zero must have an independent zero either in the same row or column; otherwise it would itself be an independent zero. Therefore each zero in a row containing no independent zeros must be crossed out vertically, as in step 112. Each independent zero must be crossed out by one line to give L

The problem of ?nding the optimal match between an N positive and an N negative sources can be solved by use of the Hungarian algorithm 1see, for example, Ref. 82. The algorithm is often discussed in terms of assigning N workers to N jobs. The cost of matching worker i to job j is ci j. The Hungarian algorithm ?nds a matching that minimizes the total cost:

C5

ooc m ,
ij ij i j

132

where mi j 5 1 if worker i is matched to job j; otherwise, mi j 5 0. Each worker must be matched to exactly one job. In this application of the Hungarian algorithm the entries in the matrix ci j 1the cost matrix2 are taken to be the distances between sources. Alternatively, the squares of these distances can be used for a minimum- sum-of-squares matching. The jobs and workers represent positive and negative sources, respectively. The following implementation of the Hungarian algorithm was used 1see, for example, Ref. 82: 112 There are N positive and N negative sources. Form an N 3 N matrix of distances between these sources. This is the cost matrix ci j. 122 Find the minimum entry in each column and subtract it from each entry in that column. 132 Find the minimum entry in each row and subtract it from each entry in that row. This is the reduced cost matrix. 142 Find the number and the position of the independent zeros in the matrix 1independent zeros are discussed below and in Subsection 2.B2. If the matrix contains N independent zeros, these represent
5102 APPLIED OPTICS @ Vol. 34, No. 23 @ 10 August 1995

crossings out in a matrix with L independent zeros. Any zeros in the same row as a vertically crossed-out independent zero cannot be crossed out horizontally and so must be crossed out vertically, as in steps 122 and 132. Further iterations of step 122 must be carried out because each vertical crossing out leads to further crossed-out independent zeros. Any 1not independent2 zeros in the same row as these newly crossed-out independent zeros will also have to be crossed out vertically. When no more crossings are made by iterations of step 122, all the remaining zeros are in the same row as an uncrossed independent zero and so can be crossed out horizontally in step 142.
B. Finding the Maximum Matching

The problem of ?nding the independent zeros in the reduced cost matrix is central to the algorithm. This is equivalent to solving the maximum-cardinality matching problem in graph theory.10 Figures 31a2 and 31b2 show a graph consisting of vertices connected by edges. In this case the vertices fall into two sets, with each edge connecting a vertex in one set 1V12 to a vertex in the other set 1V22, which represent positive and negative sources, respectively. Such a graph is known as a bipartite graph. A matching is a subset of edges with the property that no two edges are connected to the same vertex; possible matchings are shown by solid lines in Fig. 3. In a bipartite graph each vertex in V1 is connected to a vertex in V2 by at most one edge in the matching. A maximumcardinality matching is the matching that contains the most edges, subject to these constraints. Vertices not joined to an edge in the matching are known as free vertices. To relate this to the problem of ?nding the independent zeros in a matrix, we consider the zeros in a matrix represented by the edges between sets of vertices V1 and V2. It is then clear that a matching is a set of independent zeros and that ?nding a maximum-cardinality matching is equivalent to ?nding a maximal set of independent zeros. Because rows and columns of the cost matrix correspond to negative and positive sources, respectively, sources are represented by the two sets of vertices in the graphs in Fig. 3. The method of ?nding the maximum-cardinality matching, described in Ref. 10, is equivalent to ?nding the maximum ?ow in a network in the special case in which each vertex representing a positive

Fig. 3. Example of a bipartite graph. Vertices V1 and V2 represent positive and negative sources, respectively. Edges between vertices are shown by lines 1both dashed and solid2; matchings are shown by solid lines: 1a2 matching of two; 1b2 matching of three.

source is connected by an edge of unit capacity to an in?nite source; similarly, vertices representing negative sources are connected by edges of unit capacity to an in?nite sink. Edges representing connections between positive and negative sources may be taken to have in?nite capacity. The Ford–Fulkerson algorithm ?nds ?ow-augmenting paths and increases the ?ow through the network by increasing the ?ow through these paths.10 If no ?ow-augmenting path can be found, the ?ow through the network is maximal, and the algorithm terminates. A ?ow-augmenting path is equivalent in this case to a matchingaugmenting alternating path 1a path that alternately visits vertices in V1 and V2, which has the property that it starts and ends on free vertices and also has an odd number of edges2. The central part of the algorithm is to ?nd such a path. If M is the subset of edges in the graph G that form the matching between V1 and V2, then the alternating path will automatically contain edges that are alternately members or not members of M. By starting the path on an unconnected vertex, we see that the ?rst edge will not be in M and, because the path contains an odd number of edges, neither will the last edge be in M. By going along the path reversing the connections between vertices, we see that the number of connected vertices increases by two, and the matching increases by one. This is carried out to increase the matching from that in Fig. 31a2 to that in Fig. 31b2. The process is repeated for all the unconnected vertices until no matching-augmenting alternating path can be found. At this point the maximum-cardinality matching is achieved. If the maximum-cardinality matching contains fewer than N edges, where the graph represents an N 3 N matrix, then the zeros in the reduced cost matrix must be crossed out, and the matrix operations described in step 162 of the Hungarian algorithm must be carried out. However, usually only a few new zeros are created by this process, and sometimes a zero crossed in both directions is destroyed. This means that there is typically a new maximumcardinality matching that is not very different from the previous matching. It is therefore more efficient to change the matching and the graph slightly than to recalculate the graph representation of the matrix and then recalculate the maximum matching starting with zero matching: edges corresponding to zeros removed from the matrix are removed from the matching. Physically, this means that a connection between sources is discarded. If a zero is created in the matrix, the corresponding edge is added to the graph. This means that a new connection between sources is then considered. Typically the Hungarian algorithm requires more than one iteration before the maximum-cardinality matching contains N edges. At this point the matrix contains N independent zeros, and the Hungarian algorithm terminates. The positions of the independent zeros then represent the optimal matching between positive and negative sources.
10 August 1995 @ Vol. 34, No. 23 @ APPLIED OPTICS 5103

C.

Practical Implementation

The constraints of memory requirements and processing time make the direct implementation of the above algorithm impractical in its complete form. A typical noisy phase map may, for example, contain 104 sources and would therefore require a cost matrix with 108 elements, occupying at least 100 Mbytes of memory. The additional possibility of joining a source to the boundary means that the algorithm needs modi?cation in order to be applicable to the real problem. One reduces the memory requirements by considering only a small number of nearest neighbors when calculating the cost matrix. Given a source at pixel 1m, n2, a search is made for nearby sources in order of increasing distance, i.e., 1m, n 6 12 and 1m 6 1, n2, then 1m 6 1, n 6 12, followed by 1m, n 6 22 and 1m 6 2, n2, etc. Once S0 sources of opposite sign are located, the search is terminated, and the cost of joining to more distant sources is taken as in?nite. In addition, only pixels within a radius R0 of 1m, n2 are considered in order to speed up the search routine at low source density. The resulting cost matrix is sparse 1most of the elements are in?nite2 and can be stored in a more compact form. This also speeds up the matrix operations described above. Execution time for these processes becomes of order 1S0N 2 instead of order 1N 22. Values of S0 5 5 and R02 5 1000 pixels2 were used in the results presented in Section 3. The cut distributions are insensitive to the precise values of S0 and R0, although if S0 is reduced to 2 or 3, the algorithm can terminate without ?nding a minimum cost matching owing to insufficient connectivity of the graph. R0 should be at least several times the mean dipole length. Termination of the Hungarian algorithm in this reduced form is ensured provided that at no point is the minimum uncrossed entry in?nite. This is because zeros are created by subtraction of minimum uncrossed entries. If none of the minimum uncrossed entries are approximated by in?nity, then this approximation cannot affect the fact that the algorithm will terminate. If a minimum uncrossed entry is in?nite, the modi?ed Hungarian algorithm breaks down and must be rerun with larger values of S0 and@or R0. If the boundary is encountered during the search, a boundary source of opposite sign is created at the point where the boundary is detected. Typically this process results in unequal numbers of positive and negative sources, so additional ?ctitious sources are created to give equal numbers of sources of each sign. Equal numbers are needed because the Hungarian algorithm requires a square cost matrix. The cost of joining a ?ctitious source to a boundary source or to another ?ctitious source of opposite sign is taken as zero, as is the cost of joining two boundary sources. This results in unused boundary sources, being matched at zero cost. Branch cuts are placed between all the matched sources unless both sources in a given pair are boundary sources.
5104 APPLIED OPTICS @ Vol. 34, No. 23 @ 10 August 1995

Fig. 4. 1a2 Distribution of sources to be matched, showing cut lengths, with an elliptical boundary. 1b2 Optimal matching.

D.

Worked Example

A simple example illustrates the above processes. Figure 41a2 shows the positions of positive and negative sources in the wrapped phase map and the costs, i.e., the distances, or distances squared, of joining nearby sources. The sources are numbered in arbitrary order to determine the positions of their entries 1either rows or columns2 in the cost matrix. In this example, R0 5 4, so costs greater than 4 are assumed to be in?nite. S0 5 4 but does not limit the number of connections between sources to be considered in this case. Four boundary sources are created. In addition, a ?ctitious positive source 1source 72 is created to give equal numbers of positive and negative sources. This source is not shown in Fig. 4. The matrix in Fig. 51a2 shows the cost of joining sources in Fig. 41a2. The columns and rows represent positive and negative sources, respectively. S0 does not limit the possible connections between pairs of boundary sources, which can be joined at zero cost. The positive source 7 1the ?ctitious source2 may be joined only to boundary sources. The operation of subtracting the minimum entry in each column results in 2’s being subtracted from each entry in columns 1, 2, and 5, and 1’s being subtracted from columns 3 and 4. Subtracting the minimum entries in each row then gives the matrix in Fig. 51b2. The zeros in this matrix are represented in Fig. 61a2. As with Fig. 3, the vertices fall into two sets representing positive and negative sources. One pos-

Fig. 5. 1a2 Matrix representing distances between sources. Columns and rows represent 11 and 21 sources, respectively. Distances taken as in?nite are represented by dashes. 1b2 Reduced distance matrix obtained from 1a2 by subtraction of the minimum entries in each column and then in each row. Independent zeros are shown in bold; six lines are needed to cross out all the zeros. 1c2 Final matrix. The minimum uncrossed entry in 1b2, 1, is subtracted from all the uncrossed entries in 1b2 to obtain 1c2. The seven independent zeros represent the optimal matching.

zeros in the matrix are found in the same way as for the matrix of Fig. 51b2; the graph in Fig. 61b2 represents the zeros in the matrix of Fig. 51c2. The maximumcardinality matching has N 5 7, so the algorithm terminates, with the optimal matching between sources represented by the solid lines. This matching is shown in Fig. 41b2 and has a total cost of 11. The joins between unused boundary sources are not used as branch cuts acting as barriers to phase unwrapping: only joins ending on at least one real source are used.
3. Experimental

sible maximum matching is shown by the solid lines 1edges in the graph2. This is not the only possible matching containing six edges, but no matching containing seven edges is possible, so it is a maximumcardinality matching. The edges in the matching were found with the graph-theory method described above. The boldface zeros in the matrix in Fig. 51b2 correspond to the edges in the matching in the graph in Fig. 61a2. Six lines are needed to cross out all the zeros in Fig. 51b2 because there are six independent zeros. For a small matrix this can be carried out by inspection, but for larger matrices the algorithm described above to carry out step 152 of the Hungarian algorithm is used. Step 162 involves ?nding the minimum uncrossed entry, 1, and subtracting it from all the other uncrossed entries. Adding 1 to the doubly crossed entries has no effect in this example because they are already assumed to be in?nite. The matrix in Fig. 51c2 is the result of this operation. The independent

The experimental data used to test the new algorithm are the same as those described in Ref. 11, in which an in-plane speckle interferometer was used to measure in-plane rotation ?elds. Two beams from a He–Ne laser illuminated the sample at angles of 640° to the normal, giving a displacement sensitivity of ,0.5 ?m fringe21 in the horizontal in-plane direction. Images were digitized to a resolution of 512 3 512 pixels with a horizontal ?eld of view of 25 mm. Seven samples were made from aluminum sheet spray painted matte white. The resulting speckle pattern from each sample was not correlated with any of the others owing to the differing microscopic surface structure. Thus seven independent phase maps could be produced under any given conditions, permitting the accuracy and the reliability of the different algorithms to be assessed. The samples were rotated slightly about an axis perpendicular to their surfaces, and phase maps were calculated with no smoothing of the data. For a small rotation V the displacement component ux is given by

ux 5 Vy 1 u0,

142

Fig. 6. 1a2 Graph representing the zeros in the matrix in Fig. 41b2; independent zeros are represented by solid lines, other zeros by dashed lines. 1b2 Graph representing the zeros in the matrix in Fig. 41c2; solid lines represent the optimal matching.

where u0 is a constant and x and y are in-plane axes in the horizontal and vertical directions, respectively. After unwrapping, the equation for a plane 3Eq. 1424 was ?tted to the data, and the rms phase error about the best-?t solution was calculated for all the pixels in the map. This yielded a quantitative measure of the accuracy of the phase-map analysis. Three sets of experiments were done. In the ?rst set the rotation was such that 7 6 0.5 fringes spanned the ?eld of view; in the second there was a larger rotation so that there were 14 6 1 fringes; and in the third, 30 6 2 fringes. Seven phase maps were obtained from the different samples for each set. The speckle decorrelation increased with the amount of rotation; thus the number of sources grew with increasing fringe density. The phase maps were unwrapped with two new methods. Minimum-cost matching was employed to minimize the following: 1i2 the total cut length and 1ii2 the sum of squares of cut lengths. Minimizing the sum of squares has the effect of heavily penalizing long joins. The results were compared with previous methods, the most successful of which had been the modi?ed nearest-neighbor algorithm.11
10 August 1995 @ Vol. 34, No. 23 @ APPLIED OPTICS 5105

Fig. 7. Phase maps from an un?ltered speckle interferogram with 7 fringes 1left column2 and 30 fringes 1right column2 across the ?eld of view: 1a2, 1e2 wrapped maps; 1b2, 1f 2 unwrapped without branch cuts; 1c2, 1g2 unwrapped with branch cuts placed by modi?ed nearest-neighbor method; 1d2, 1h2 unwrapped with minimum-costmatching branch cuts. The linear ramp is subtracted from 1c2, 1d2, 1g2, and 1h2 to show the absolute value of the phase error. White and black represent zero and 4p phase errors, respectively.

Some typical wrapped phase maps for the lowest and highest fringe densities are shown in Figs. 71a2 and 71e2. These are raw phase maps, with no smoothing of the data. The conventional method of unwrapping that uses no branch cuts is unsuccessful 3see Figs. 71b2 and 71f 24. The unwrapped map does not look like a plane, and the rms error is tens of times

larger than that obtained from the other methods. Figures 71c2 and 71g2 were unwrapped by the modi?ed nearest-neighbor method, and Figs. 71d2 and 71h2 were unwrapped by the minimum-cost-matching algorithm. In all four cases the image shows the residuals between the unwrapped map and the best-?t plane in order to enhance the contrast and to highlight any unwrapping errors. Minimum-cost matching was the only method to work reliably with the 30-fringe data 1,0.1 sources pixel212. The other methods all failed by producing long, physically unreasonable branch cuts and corresponding discontinuities in the unwrapped phase map. Cut distributions from unwrapping the top-right corner of Fig. 71e2 1one third of both the horizontal and the vertical ?elds of view2 by three methods are shown in Fig. 8, and rms phase errors are compared graphically in Fig. 9. The error bars represent the standard deviation in rms phase error among the values in each set. The nearest-neighbor method fails at source densities of 0.05 sources pixel21, whereas the modi?ed nearest-neighbor method fails at 0.1 sources pixel21. The minimum-cost-matching method gives the best performance at high source densities 1rms phase error of 1.2 rad, compared with 1.8 rad for the modi?ed nearest-neighbor method at 0.1 sources pixel212, and it also gives slightly smaller errors at low densities. At high source density, minimizing the sum of squares of the cut lengths gives a small 13%2 improvement over minimizing the total cut length, but there is also a reduction in the standard deviation of the rms phase errors, showing that the minimum-squares method is slightly more consistent. The rms phase error of 1.2 rad is 0.6% of the full phase change over the image 1190 rad2. The algorithm was executed on a Sun SPARCstation 2 1benchmark ratings of 28.5 Mips and 4.2 M?ops2. Typical elapsed clock times for the branchcut placement routines were, on a map containing 22,000 sources, ,60 s for the nearest-neighbor and the modi?ed nearest-neighbor algorithms, 400 s for the minimum-cost-matching algorithm when the total cut length is minimized, and 700 s when the sum of the cut length squared is minimized. Computation times for both modi?ed nearest-neighbor and minimum total-cut-length methods increase approximately as N 2.0, where N is the number of sources.

Fig. 8. Cut distributions for a phase map consisting of the top one ninth of Fig. 71e2 placed by 1a2 the nearest-neighbor method, 1b2 the modi?ed nearest-neighbor method, and 1c2 the minimum-cost-matching algorithm.
5106 APPLIED OPTICS @ Vol. 34, No. 23 @ 10 August 1995

cally unreasonable 2p discontinuities in the unwrapped map 3Fig. 101b24. The minimum-cost-matching algorithm has so far been used successfully on ,102 phase maps and has placed a total of 105–106 branch cuts with no discernible unwrapping errors, provided R0 and S0 are sufficiently large.
4. Conclusions

Fig. 9. Rms phase error obtained by unwrapping with four branch-cut placement algorithms on speckle-interferometry data with three different discontinuity source densities.

The time for minimum-cost matching increased as N 2.6 when the sum of the cut length squared was minimized, because a larger number of iterations of the Hungarian algorithm were needed. The source densities used above are higher than would be found in many practical applications, where low-pass ?ltering can be used to reduce the source density in speckle interferograms. However, in some situations a local region of low SNR can cause global failure of all but the minimum-cost-matching method. The example in Fig. 10 shows a disk of a ?lled polymer 1diameter of 10 mm2 being loaded vertically in compression. The wrapped phase map 3Fig. 101a24 represents the horizontal in-plane displacement when compressed, with a sensitivity of 0.433 ?m per fringe, and was measured by phase-stepping speckle interferometry. Many of the discontinuity sources were removed by prior smoothing of the correlation fringe patterns with a 3 3 3 ?lter. However, speckle decorrelation associated with the formation of a tear near the center of the disk causes failure of all methods except the minimum-cost-matching algorithm. Figures 101b2 and 101c2 show the result of unwrapping Fig. 101a2 by the modi?ed nearestneighbor method and the minimum-cost-matching algorithm, respectively. Long cuts appear as physi-

A new method 1the minimum-cost-matching method2 for placing branch cuts for phase unwrapping of noisy phase maps has been developed. The new method works reliably at source densities of 0.1 sources pixel21, at which the previous best 1modi?ed nearestneighbor2 method failed. The minimum-cost-matching method uses a graph-theory approach to ?nd the globally optimal placement of branch cuts, minimizing the total cut length in time of order 1N 22. This makes a global optimization practical: a direct evaluation of all possible matchings would require time of order N !, and a simulated-annealing approach to global optimization was found to be very slow.12 The minimum-cost-matching method requires only minutes of computation time on a relatively inexpensive workstation even for the noisiest of phase maps, containing upwards of 104 sources. It is therefore recommended as the best of the branch-cut placement algorithms. The method is directly applicable to noisy phase maps in which the underlying phase distributions are continuous across the ?eld of view. When real phase discontinuities arise 1owing to steps, cracks, specimen edges, etc.2, these must be represented by boundaries; otherwise, isolated discontinuity sources can occur that result in the global failure of the branch-cut method. Simply connected regions can then be unwrapped reliably by direct application of the proposed algorithm. Multiply connected regions 1i.e., regions of the image containing separate boundaries within boundaries2 normally can also be unwrapped correctly, but an unacceptable cut distribution could be produced because there is no mechanism to constrain the sum of the sources joined to a boundary to be zero. However, this problem could in principle be overcome simply by linking of the internal boundaries with cuts to convert the multiply connected

Fig. 10. 1a2 Wrapped phase map for a disk undergoing compression between two anvils; phase maps from 1a2, unwrapped by 1b2 the modi?ed nearest-neighbor method and 1c2 the minimum-cost-matching method. The magnitude of the phase is shown 1negative phases displayed as positive2 to enhance contrast.
10 August 1995 @ Vol. 34, No. 23 @ APPLIED OPTICS 5107

region into a simply connected region, as in the analysis of functions of a complex variable. The main remaining obstacle to a robust, fully automated phase-unwrapping algorithm is the problem of setting up the boundaries automatically. This is not, in general, possible, given just a single phase map, owing to the problem of uniqueness: many different unwrapped maps, with true phase discontinuities in different places, could give rise to the same wrapped phase map. Future progress in this area is therefore likely to depend on incorporating additional information, such as phase maps recorded at different times or with different interferometric sensitivities, or on independent measurements of the specimen geometry. We are grateful to H. T. Goldrein for help with the development of the speckle interferometer. J. M. Huntley is pleased to acknowledge a research fellowship from the Royal Society, and S. R. E. Turner was supported by a research studentship from the Engineering and Physical Sciences Research Council.
References
1. G. T. Reid, ‘‘Automatic fringe pattern analysis: Opt. Lasers Eng. 7, 37–68 119862. a review,’’

2. D. W. Robinson, ‘‘Phase unwrapping methods,’’ Interferogram Analysis, D. W. Robinson and G. T. Reid, eds. 1Institute of Physics, Bristol, UK, 19932, pp. 194–229. 3. R. M. Goldstein, H. A. Zebker, and C. L. Werner, ‘‘Satellite radar interferometry: two-dimensional phase unwrapping,’’ Radio Sci. 23, 713–720 119882. 4. J. M. Huntley, ‘‘Noise-immune phase unwrapping algorithm,’’ Appl. Opt. 28, 3268–3270 119892. 5. L. D. Barr, V. Coude ? du Foresto, J. Fox, G. A. Poczulp, J. Richardson, C. Roddier, and F. Roddier, ‘‘Large-mirror testing facility at the National Optical Astronomy Observatories,’’ Opt. Eng. 30, 1405–1414 119912. 6. J. M. Huntley, ‘‘New methods for unwrapping noisy phase maps,’’ presented at the SPIE International Conference on Interferometry, Warsaw, Poland, 16–20 May, 1994. 7. J. M. Huntley and J. R. Buckland, ‘‘Characterization of sources of 2p phase discontinuity in speckle interferograms,’’ J. Opt. Soc. Am. A 1to be published2. 8. D. R. Wilson, ‘‘Minimax theorems in graph theory,’’ Selected Topics in Graph Theory, L. W. Beineke and R. J. Wilson, eds. 1Academic, New York, 19782, pp. 250–252. 9. M. S. Bazaraa, J. J. Jarvis, and H. D. Sherali, Linear Programming and Network Flows, 2nd ed. 1Wiley, New York, 19902. 10. A. Gibbons, Algorithmic Graph Theory, 1Cambridge U. Press, Cambridge, UK, 19852. 11. R. Cusack, J. M. Huntley, and H. T. Goldrein, ‘‘Improved noise-immune phase unwrapping algorithm,’’ Appl. Opt. 34, 781–789 119952.

5108

APPLIED OPTICS @ Vol. 34, No. 23 @ 10 August 1995


相关文章:
更多相关标签: