Abstract

Two-dimensional (2D) phase unwrapping continues to find applications in a wide variety of scientific and engineering areas including optical and microwave interferometry, adaptive optics, compensated imaging, and synthetic-aperture-radar phase correction, and image processing. We have developed a robust method (not based on any path-following scheme) for unwrapping 2D phase principal values (in a least-squares sense) by using fast cosine transforms. If the 2D phase values are associated with a 2D weighting, the fast transforms can still be used in iterative methods for solving the weighted unwrapping problem. Weighted unwrapping can be used to isolate inconsistent regions (i.e., phase shear) in an elegant fashion.

© 1994 Optical Society of America

Full Article  |  PDF Article

References

  • View by:
  • |
  • |
  • |

  1. D. L. Fried, “Least-squares fitting a wave-front distortion estimate to an array of phase-difference measurements,”J. Opt. Soc. Am. 67, 370–375 (1977).
    [CrossRef]
  2. R. H. Hudgin, “Wave-front reconstruction for compensated imaging,”J. Opt. Soc. Am. 67, 375–378 (1977).
    [CrossRef]
  3. R. J. Noll, “Phase estimates from slope-type wave-front sensors,”J. Opt. Soc. Am. 68, 139–140 (1978).
    [CrossRef]
  4. B. R. Hunt, “Matrix formulation of the reconstruction of phase values from phase differences,”J. Opt. Soc. Am. 69, 393–399 (1979).
    [CrossRef]
  5. H. Takajo, T. Takahashi, “Least-squares phase estimation from phase differences,” J. Opt. Soc. Am. A 5, 416–425 (1988).
    [CrossRef]
  6. H. T. Takajo, T. Takahashi, “Noniterative method for obtaining the exact solution for the normal equation in least-squares phase estimation from the phase difference,” J. Opt. Soc. Am. A 5, 1818–1827 (1988).
    [CrossRef]
  7. D. C. Ghiglia, G. A. Mastin, “Two-dimensional phase correction of synthetic-aperture-radar imagery,” Opt. Lett. 14, 1104–1106 (1989).
    [CrossRef] [PubMed]
  8. D. C. Ghiglia, L. A. Romero, “Direct phase estimation from phase differences using fast elliptic partial differential equation solvers,” Opt. Lett. 14, 1107–1109 (1989).
    [CrossRef] [PubMed]
  9. B. L. Busbee, G. H. Gollub, C. W. Nielson, “On direct methods for solving Poisson’s equations,” SIAM J. Numer. Anal. 7, 627–656 (1970).
    [CrossRef]
  10. D. C. Ghiglia, G. A. Mastin, L. A. Romero, “Cellular automata method for phase unwrapping,” J. Opt. Soc. Am. A 4, 267–280 (1987).
    [CrossRef]
  11. R. M. Goldstein, H. A. Zebker, C. L. Werner, “Satellite radar interferometry: Two-dimensional phase unwrapping,” Radio Sci. 23, 713–720 (1988).
    [CrossRef]
  12. N. H. Ching, D. Rosenfeld, M. Braun, “Two-dimensional phase unwrapping using a minimum spanning tree algorithm,”IEEE Trans. Image Process. 1, 355–365 (1992).
    [CrossRef] [PubMed]
  13. J. S. Lim, “The discrete cosine transform,” in Two-Dimensional Signal and Image Processing (Prentice-Hall, Englewood Cliffs, N.J., 1990), pp. 148–157.
  14. G. H. Gollub, C. F. Van Loan, “Iterative methods for linear systems,” in Matrix Computations, 2nd ed. (Johns Hopkins U. Press, Baltimore, Md., 1990), pp. 516–538.
  15. cosft must be used in the forward mode in all cases. In addition, the technique mentioned on p. 649 of Ref. 16 for solving Poisson’s equation with Neumann boundary conditions by using cosft alone is not correct and does not work. cosft does not specify the boundary conditions exactly [as defined by Eq. (14)] and therefore does not solve the least-squares problem.
  16. W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling, Numerical Recipes: The Art of Scientific Computing (Cambridge U. Press, Cambridge, 1986).
  17. G. A. Mastin, D. C. Ghiglia, “A research-oriented spotlight synthetic aperture radar polar reformatter,” Publ. SAND90-1793 (Sandia National Laboratories, Albuquerque, N.M., 1990), pp. 23–27.

1992 (1)

N. H. Ching, D. Rosenfeld, M. Braun, “Two-dimensional phase unwrapping using a minimum spanning tree algorithm,”IEEE Trans. Image Process. 1, 355–365 (1992).
[CrossRef] [PubMed]

1989 (2)

1988 (3)

1987 (1)

1979 (1)

1978 (1)

1977 (2)

1970 (1)

B. L. Busbee, G. H. Gollub, C. W. Nielson, “On direct methods for solving Poisson’s equations,” SIAM J. Numer. Anal. 7, 627–656 (1970).
[CrossRef]

Braun, M.

N. H. Ching, D. Rosenfeld, M. Braun, “Two-dimensional phase unwrapping using a minimum spanning tree algorithm,”IEEE Trans. Image Process. 1, 355–365 (1992).
[CrossRef] [PubMed]

Busbee, B. L.

B. L. Busbee, G. H. Gollub, C. W. Nielson, “On direct methods for solving Poisson’s equations,” SIAM J. Numer. Anal. 7, 627–656 (1970).
[CrossRef]

Ching, N. H.

N. H. Ching, D. Rosenfeld, M. Braun, “Two-dimensional phase unwrapping using a minimum spanning tree algorithm,”IEEE Trans. Image Process. 1, 355–365 (1992).
[CrossRef] [PubMed]

Flannery, B. P.

W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling, Numerical Recipes: The Art of Scientific Computing (Cambridge U. Press, Cambridge, 1986).

Fried, D. L.

Ghiglia, D. C.

Goldstein, R. M.

R. M. Goldstein, H. A. Zebker, C. L. Werner, “Satellite radar interferometry: Two-dimensional phase unwrapping,” Radio Sci. 23, 713–720 (1988).
[CrossRef]

Gollub, G. H.

B. L. Busbee, G. H. Gollub, C. W. Nielson, “On direct methods for solving Poisson’s equations,” SIAM J. Numer. Anal. 7, 627–656 (1970).
[CrossRef]

G. H. Gollub, C. F. Van Loan, “Iterative methods for linear systems,” in Matrix Computations, 2nd ed. (Johns Hopkins U. Press, Baltimore, Md., 1990), pp. 516–538.

Hudgin, R. H.

Hunt, B. R.

Lim, J. S.

J. S. Lim, “The discrete cosine transform,” in Two-Dimensional Signal and Image Processing (Prentice-Hall, Englewood Cliffs, N.J., 1990), pp. 148–157.

Mastin, G. A.

Nielson, C. W.

B. L. Busbee, G. H. Gollub, C. W. Nielson, “On direct methods for solving Poisson’s equations,” SIAM J. Numer. Anal. 7, 627–656 (1970).
[CrossRef]

Noll, R. J.

Press, W. H.

W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling, Numerical Recipes: The Art of Scientific Computing (Cambridge U. Press, Cambridge, 1986).

Romero, L. A.

Rosenfeld, D.

N. H. Ching, D. Rosenfeld, M. Braun, “Two-dimensional phase unwrapping using a minimum spanning tree algorithm,”IEEE Trans. Image Process. 1, 355–365 (1992).
[CrossRef] [PubMed]

Takahashi, T.

Takajo, H.

Takajo, H. T.

Teukolsky, S. A.

W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling, Numerical Recipes: The Art of Scientific Computing (Cambridge U. Press, Cambridge, 1986).

Van Loan, C. F.

G. H. Gollub, C. F. Van Loan, “Iterative methods for linear systems,” in Matrix Computations, 2nd ed. (Johns Hopkins U. Press, Baltimore, Md., 1990), pp. 516–538.

Vetterling, W. T.

W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling, Numerical Recipes: The Art of Scientific Computing (Cambridge U. Press, Cambridge, 1986).

Werner, C. L.

R. M. Goldstein, H. A. Zebker, C. L. Werner, “Satellite radar interferometry: Two-dimensional phase unwrapping,” Radio Sci. 23, 713–720 (1988).
[CrossRef]

Zebker, H. A.

R. M. Goldstein, H. A. Zebker, C. L. Werner, “Satellite radar interferometry: Two-dimensional phase unwrapping,” Radio Sci. 23, 713–720 (1988).
[CrossRef]

IEEE Trans. Image Process. (1)

N. H. Ching, D. Rosenfeld, M. Braun, “Two-dimensional phase unwrapping using a minimum spanning tree algorithm,”IEEE Trans. Image Process. 1, 355–365 (1992).
[CrossRef] [PubMed]

J. Opt. Soc. Am. (4)

J. Opt. Soc. Am. A (3)

Opt. Lett. (2)

Radio Sci. (1)

R. M. Goldstein, H. A. Zebker, C. L. Werner, “Satellite radar interferometry: Two-dimensional phase unwrapping,” Radio Sci. 23, 713–720 (1988).
[CrossRef]

SIAM J. Numer. Anal. (1)

B. L. Busbee, G. H. Gollub, C. W. Nielson, “On direct methods for solving Poisson’s equations,” SIAM J. Numer. Anal. 7, 627–656 (1970).
[CrossRef]

Other (5)

J. S. Lim, “The discrete cosine transform,” in Two-Dimensional Signal and Image Processing (Prentice-Hall, Englewood Cliffs, N.J., 1990), pp. 148–157.

G. H. Gollub, C. F. Van Loan, “Iterative methods for linear systems,” in Matrix Computations, 2nd ed. (Johns Hopkins U. Press, Baltimore, Md., 1990), pp. 516–538.

cosft must be used in the forward mode in all cases. In addition, the technique mentioned on p. 649 of Ref. 16 for solving Poisson’s equation with Neumann boundary conditions by using cosft alone is not correct and does not work. cosft does not specify the boundary conditions exactly [as defined by Eq. (14)] and therefore does not solve the least-squares problem.

W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling, Numerical Recipes: The Art of Scientific Computing (Cambridge U. Press, Cambridge, 1986).

G. A. Mastin, D. C. Ghiglia, “A research-oriented spotlight synthetic aperture radar polar reformatter,” Publ. SAND90-1793 (Sandia National Laboratories, Albuquerque, N.M., 1990), pp. 23–27.

Cited By

OSA participates in CrossRef's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (12)

Fig. 1
Fig. 1

Wrapped values of a 512 × 512-pixel 2D phase plane (scaled for display).

Fig. 2
Fig. 2

Unwrapped phase from values depicted in Fig. 1.

Fig. 3
Fig. 3

Rewrapped values from Fig. 2 for direct comparison with Fig. 1. In this case the least-squares unwrapping of totally consistent data yielded perfect results.

Fig. 4
Fig. 4

Wrapped 2D planar phase values with a small rectangular region containing uniformly distributed noise. The phase is inconsistent within the rectangular region and around its boundary.

Fig. 5
Fig. 5

Rewrapped values from algorithm 1 applied to Fig. 4. Inconsistent data influences the unwrapped results locally, with diminishing effects further away.

Fig. 6
Fig. 6

Unwrapped values of Fig. 4 from algorithm 1 showing minimal influence of noisy region on the final result.

Fig. 7
Fig. 7

Two-dimensional weighting array to accompany the wrapped phase data depicted in Fig. 4. Black corresponds to zero-valued weights, and white corresponds to unity weight. The dark border defines the array boundary.

Fig. 8
Fig. 8

Sequence of Picard iteration weighted phase unwrapping (algorithm 2) applied to Fig. 4 with the corresponding weighting array of Fig. 7. (a) One iteration, (b) two iterations, (c) three iterations, (d) ten iterations. Note correct solution in region corresponding to unity weights and continuity across region corresponding to zero-valued weights. We show rewrapped phases to highlight subtle features.

Fig. 9
Fig. 9

Wrapped noise-free phase values depicting a phase surface with a shear along the horizontal line midway between the top and the bottom of the figure.

Fig. 10
Fig. 10

Unweighted 2D least-squares unwrap of phases that depict a phase shear. (a) Unwrapped result, (b) rewrapped values shown to illustrate solution subtleties and least-squares seaming across the horizontal shear line.

Fig. 11
Fig. 11

Weighted unwrap (algorithm 3) of phase data with shear. (a) Rewrapped output after 10 iterations, (b) rewrapped output after 20 iterations, and (c) unwrapped output after 20 iterations. Note that top and bottom halves unwrapped independently, with no influence across the shear line.

Fig. 12
Fig. 12

Some rewrapped results of algorithm 2 applied to Fig. 9 for comparison with the convergence rate of algorithm 3. (a) 10 iterations, (b) 20 iterations, (c) 30 iterations, and (d) 100 iterations. Final convergence was achieved at approximately 500 iterations, in contrast to the convergence at 20 iterations for algorithm 3.

Equations (39)

Equations on this page are rendered with MathJax. Learn more.

ψ i , j = ϕ i , j + 2 π k , k an integer , - π < ψ i , j π , i = 0 M - 1 , j = 0 N - 1.
Δ i , j x = W { ψ i + 1 , j - ψ i , j } i = 0 M - 2 , j = 0 N - 1 , Δ i , j x = 0 , otherwise ;
Δ i , j y = W { ψ i , j + 1 - ψ i , j } , i = 0 M - 1 , j = 0 N - 2 , Δ i , j y = 0 , otherwise ,
i = 0 M - 2 j = 0 N - 1 ( ϕ i + 1 , j - ϕ i , j - Δ i , j x ) 2 + i = 0 M - 1 j = 0 N - 2 ( ϕ i , j + 1 - ϕ i , j - Δ i , j y ) 2
ϕ i + 1 , j + ϕ i - 1 , j + ϕ i , j + 1 + ϕ i , j - 1 - 4 ϕ i , j = Δ i , j x - Δ i - 1 , j x + Δ i , j y - Δ i , j - 1 y .
( ϕ i + 1 , j - 2 ϕ i , j + ϕ i - 1 , j ) + ( ϕ i , j + 1 - 2 ϕ i , j + ϕ i , j - 1 ) = ρ i , j ,
ρ i , j = ( Δ i , j x - Δ i - 1 , j x ) + ( Δ i , j y - Δ i , j - 1 y ) .
2 x 2 ϕ ( x , y ) + 2 y 2 ϕ ( x , y ) = ρ ( x , y ) .
Δ - 1 , j x = 0 Δ M - 1 , j x = 0 , j = 0 N - 1 ,
Δ i , - 1 y = 0 , Δ i , N - 1 y = 0 , i = 0 M - 1.
C m , n = { i = 0 M - 1 j = 0 N - 1 4 x i , j cos [ π 2 M m ( 2 i + 1 ) ] cos [ π 2 N n ( 2 j + 1 ) ] 0 m M - 1 ;             0 n N - 1 0 otherwise .
x i , j = { 1 M N m = 0 M - 1 n = 0 N - 1 w 1 ( m ) w 2 ( n ) C m , n cos [ π 2 M m ( 2 i + 1 ) ] cos [ π 2 N n ( 2 j + 1 ) ] 0 i M - 1 ,             0 j N - 1 0 otherwise , w 1 ( m ) = 1 / 2 , m = 0 , w 1 ( m ) = 1 , 1 m M - 1 , w 2 ( m ) = 1 / 2 , n = 0 , w 2 ( m ) = 1 , 1 n N - 1.
ϕ i , j = 1 M N m = 0 M - 1 n = 0 N - 1 w 1 ( m ) w 2 ( n ) ϕ ^ m , n × cos [ π 2 M m ( 2 i + 1 ) ] cos [ π 2 N n ( 2 j + 1 ) ]
ϕ ^ i , j = ρ ^ i , j 2 ( cos π i M + cos π j N - 2 ) .
ϕ 0 , j - ϕ - 1 , j = 0 , ϕ M , j - ϕ M - 1 , j = 0 , j = 0 N - 1 ϕ i , 0 - ϕ i , - 1 = 0 , ϕ i , N - ϕ i , N - 1 = 0 , i = 0 M - 1.
Ax = b ,
A T Ax = A T b ,
P ϕ = ρ ,
WAx = Wb .
A T W T WAx = A T W T Wb ,
Qx = A T b ¯ .
c = A T b ¯ ,
Q ϕ = c .
Q = P + D .
( P + D ) ϕ = c
P ϕ = c - D ϕ .
P ϕ k + 1 = c - D ϕ k ,
D ϕ k = ( Q - P ) ϕ k ,
β k = r k - 1 T z k - 1 / r k - 2 T z k - 2 , p k = z k - 1 + β k p k - 1 .
α k = r k - 1 T z k - 1 / p k T Qp k , ϕ k = ϕ k - 1 + α k p k , r k = r k - 1 - α k Qp k .
C m = { i = 0 M - 1 2 x i cos [ π 2 M m ( 2 i + 1 ) ] 0 m M - 1 0 otherwise ,
x i = { 1 M m = 0 M - 1 w 1 ( m ) C m [ cos π 2 M m ( 2 i + 1 ) ] 0 i M - 1 0 otherwise , w 1 ( m ) = 1 / 2 , m = 0 , w 1 ( m ) = 1 , 1 m M - 1.
C m = 2 cos π m 2 M i = 0 M - 1 x i cos π m M i - 2 sin π m 2 M i = 0 M - 1 x i sin π m M i ,
x i = 1 M m = 0 M - 1 C ^ 1 m cos π m M i - 1 M i = 0 M - 1 C ^ 2 m sin π m M i ,
C ^ 1 m = w 1 ( m ) C m cos π m 2 M ,
C ^ 2 m = w 1 ( m ) C m sin π m 2 M .
c i , j = min ( w i + 1 , j 2 , w i , j 2 ) Δ i , j x - min ( w i , j 2 , w i - 1 , j 2 ) Δ i - 1 , j x + min ( w i , j + 1 2 , w i , j 2 ) Δ i , j y - min ( w i , j 2 , w i , j - 1 2 ) Δ i , j - 1 y .
ρ i , j = c i , j - [ ( w x 1 2 - 1 ) ( ϕ i + 1 , j - ϕ i , j ) - ( w x 2 2 - 1 ) ( ϕ i , j - ϕ i - 1 , j ) + ( w y 1 2 - 1 ) ( ϕ i , j + 1 - ϕ i , j ) - ( w y 2 2 - 1 ) ( ϕ i , j - ϕ i , j - 1 ) ] ,
w x 1 2 = min ( w i + 1 , j 2 , w i , j 2 ) , w x 2 2 = min ( w i , j 2 , w i - 1 , j 2 ) , w y 1 2 = min ( w i , j + 1 2 , w i , j 2 ) , w y 2 2 = min ( w i , j 2 , w i , j - 1 2 ) ,

Metrics