# A Novel Phase-Based Approach to Tear Film Surface Quality Assessment Using Lateral Shearing Interferometry

Piotr Szyperski, D. Robert Iskander

<figure><img src="image_1.png" /></figure>

## To cite this version:

Piotr Szyperski, D. Robert Iskander. A Novel Phase-Based Approach to Tear Film Surface Quality Assessment Using Lateral Shearing Interferometry. 14th Computer Information Systems and Industrial Management (CISIM), Sep 2015, Warsaw, Poland. pp.435-447, 10.1007/978-3-319-24369-6_36 hal-01444485

HAL Id: hal-01444485

https://inria.hal.science/hal-01444485

Submitted on 24 Jan 2017

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A novel phase-based approach
to tear film surface quality assessment
using lateral shearing interferometry

Piotr Szyperski* and D. Robert Iskander

Wrocław University of Technology, Department of Computer Engineering,
Wybrzeże Wyspiańskiego 27, 50-370 Wrocław, Poland,

{piotr.szyperski,robert.iskander}@pwr.edu.pl,

http://www.tearfilm.pwr.wroc.pl/index.php/en/research,

http://dri.pwr.edu.pl/research/

**Abstract.** Lateral shearing interferometry (LSI) can be used for assessing the properties of tear film. In particular, it has the capability of acquiring dynamic variations in tear film surface quality (TFSQ) during an interblink interval in an *in-vivo* fashion. The purpose of this study was to assess the suit-ability of LSI based two-dimensional (2-D) phase estimation procedures in the analysis of tear film dynamics. The paper discusses the main difficulty in 2-D phase estimation – the problem of phase unwrapping, proposes a modi-fication of one of the popular phase unwrapping algorithms, and suggests a set of phase parameters that could be exploited as LSI-based TFSQ descriptors.

**Keywords:** *in-vivo* interferometry, Takeda's algorithm, two-dimensional phase unwrapping, Goldstein's algorithm.

# 1 Introduction

*Background.* Lateral shearing interferometry (LSI) is one of the most sensitive me-thods for measuring tear film surface quality (TFSQ) *in-vivo* in human eyes [1]. Current analysis of LSI-based interferograms essentially includes parameters derived based on the analysis of the first harmonic of the amplitude Fourier spectrum [2]. It was of interest to ascertain whether the phase contained in the spectrum can be used to extract additional information from the interferogram that can be further utilized for the assessment of TFSQ. A classical method for phase estimation in interferometry and profilometry is based on the Takeda approach [3]. The main problem in LSI is that the spatial phase requires two-dimensional (2-D) unwrapping, for which the algorithms are not trivial [4] and computationally very intensive. Furthermore, the results of 2-D phase unwrapping are highly dependent on the particular algorithm used. The aim of this work was to assess the applicability of 2-D phase estimation for LSI and design a set of phase parameters for the characterization of the dynamics in TFSQ analysis both on the eye's cornea and contact lens surface.

<footnote>* Corresponding author.</footnote>
## 1.1 Hardware Setup of Experiments

The LSI hardware setup under consideration comprises single set optical wedge with predetermined shear angle. Conversely, Twyman-Green utilizes the reference mirror that allows measuring wavefront deformations in a direct manner [5–7]. The drawback of this setup reducing its applicability in TFSQ assessment is the fact that it is more sensitive to even slight vibrations of the test surface (e. g. saccadic movements of the eye). These are inevitable during *in-vivo* experiments whereas LSI experimental setup is capable of mitigating their impact on the interference pattern [2].

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** LSI scheme: M – mirrors, BS – beam splitter, OW – optical wedge</figcaption></figure>

Figure 1 depicts a typical LSI configuration employed for corneal and contact lens tear film surface topography data acquisition while Fig. 2 shows few examples of output interferograms.

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** typical cases of tear film interferograms showing pre-corneal tear film in good condition (a), slightly deteriorated (b) and pre-lens tear film (c)</figcaption></figure>
Until now the analysis of the LSI images concerned only the amplitude information and little attention was paid to the phase [8]. The purpose of this work was to examine the applicability of phase-based methods in the analysis of tear film.

## 2 Methods

The true wavefront phase cannot be estimated directly from LSI images [9], although some methods of evaluating an unknown wavefronts have been developed, e. g. Saunders method or Rimmer-Wyant method with Okuda improvements. The main limitation of these methods is that the wavefront evaluation is confined only to points separated by a distance of *S* that is the lateral shear in the sagittal (0*x*) direction [7, Sec. 1.5.2, 1.5.3].

### 2.1 Indirect Retrieval of Phase Information

There is some visual similarity between interferometric images in the single direction LSI and those of *profilometry*. Hence, it was of interest to ascertain whether methods developed for phase estimation in profilometry would be applicable also to LSI. The basic principle of phase estimation in profilometry has been outlined by Takeda [3, 10] in which both 1-D and 2-D phase estimation algorithms were considered. In profilometry, such algorithms are usually followed by a phase unwrapping procedure which reveals the true surface of the wavefront exhibiting the shape of the scanned object. In case of 3-D objects, an algorithm for 2-D phase unwrapping has been proposed by Goldstein's *et al.* [11] for application in satellite radar interferometry to measure ground surface deformation or terrain elevation. Goldstein's algorithm is considered as a good trade-off between trivial direct phase unwrapping methods and more sophisticated but also computationally very expensive procedures [4, 12]. Furthermore, for which there are few ready implementation routines.

**The Takeda's Phase Reconstruction Algorithm.** This method was proposed by M. Takeda *et al.* in 1981 to enable automatic discrimination between elevation and depression of the object or wavefront form [3]. The chain of the operations constituting the routine is depicted by Fig. 3 and discussed below in greater detail.

<figure><img src="image_5.png" /><figcaption>Fig. 3. image processing chain in Takeda’s algorithm</figcaption></figure>
Assumptions on fringe pattern periodicity and phase disruptions. When a periodic fringe pattern is considered, it can be represented in the form of:

$$g(x,y)=a(x,y)+b(x,y)\cos[2\pi f_0 x+\phi(x,y)], \quad (1)$$

where $\phi(x,y)$ contains the desired phase information and $a(x,y)$ and $b(x,y)$ are unwanted irradiance variations that, together with $\phi(x,y)$, in most cases present much slower spatial frequency than the variation introduced by $f_0$. Equation (1) can be rewritten equivalently as:

$$g(x,y)=a(x,y)+c(x,y)\exp(2\pi i f_0 x)+c^*(x,y)\exp(-2\pi i f_0 x), \quad (2)$$

where

$$c(x,y)=\frac{1}{2}b(x,y)\exp[i\phi(x,y)], \quad (3)$$

with $*$ denoting a complex conjugate.

Carrier sideband detection and extraction in spectral domain. During the analysis of 2-D interferogram, the following procedure can be performed in either of the two dimensions as well as in both. Since the operation in 2-D do not differ significantly from 1-D case, the latter is described to maintain simplicity of the notation. The Fast Fourier Transformation (FFT) algorithm can be applied to the fringe pattern along the axis of abscissæ. The result of transforming Eq. (2) into spatial frequency domain is:

$$G(f,y)=A(f,y)+C(f-f_0,y)+C^*(f+f_0,y), \quad (4)$$

with capital letters denoting the Fourier spectra (and $C$ used for the carrier side bands) and $f$ being the spatial frequency in $0x$ direction.

In a typical case of an interferogram with low distortion grade, the sidebands Fourier spectra from Eq. (4) would be well-separated from the carrier frequency $f_0$. This is ensured as long as the variations of $a(x,y)$, $b(x,y)$ and $\phi(x,y)$ are sufficiently slower than the spatial frequency $f_0$. Either of the two sidebands spectra on the carrier can be used for phase information retrieval, and the commonly used is $C(f-f_0,y)$.

A few simple operations are applied on the periodogram to mask and extract the carrier sideband spectrum. Then it is translated by $-f_0$ on the frequency axis toward its origin to obtain $C(f,y)$. Note that the unwanted background variation $a(x,y)$ representing has been filtered out in this stage. As $A(f,y)$ represents the value of the 0th order harmonic component in Fourier spectrum, this reduction can be further improved by shifting the interferogram brightness values during the preprocessing to reduce its mean to have zero value.

Fringe phase reconstruction. The Inverse FFT (IFFT) algorithm is used to compute $c(x,y)$, the inverse Fourier transform of $C(f,y)$ with respect to $f$. It is defined by Eq. (3) and a complex natural logarithm of it can be computed:

$$\ln[c(x,y)]=\ln[\frac{1}{2}b(x,y)]+i\phi(x,y). \quad (5)$$

The phase $\phi(x,y)$ constitutes the imaginary part of the result, as Eq. (5) shows. Furthermore, it is completely separated from the real part that comprises the unwanted amplitude variation $b(x,y)$. The phase obtained in this manner is in fact
the angular coordinate of $C(f,y)$ represented in polar form. Thus, it is indeterminate to a factor of $2\pi$, which is often referred to as *wrapped phase*. In most cases, a computer-generated function subroutine gives a principal value ranging from $-\pi$ to $\pi$. This can lead to presence of discontinuities in reconstructed phase that sometimes need to be eliminated to allow further analysis of the absolute phase value in a process commonly referred to as *phase unwrapping* that is described further.

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Gaussian phase disruption to be reconstructed by Takeda's algorithm: (a) unwrapped, used for phase modulation of fringes, and (b) wrapped mod $2\pi$ for mimicking the error-less output (disruption scale increases from left to right)</figcaption></figure>

*Advantages and Limitations of Takeda's Algorithm.* The method described above is more accurate and sensitive than e. g. moiré topography or conventional fringe-contour interferometry and does not require any moving components in experimental setup nor involve interpolation of the data in the regions between contour fringes [3, 10]. Nevertheless, under specific circumstances, its results can be inaccurate. This was investigated on synthetic data involving Gaussian phase disruptions of maximum intensities of $\phi_{max}(x,y)$ varying from low to high as well as low, intermediate and high fringe frequencies $f_0$.

The phase disruption surface used for investigation of Takeda's algorithm implementation has maximum value in the center and zeros in the corners, q. v. Fig. 4. The wrapped phase surfaces in Fig. 4b were generated for comparison with actual results. Figures 5 to 8 show (clockwise): Takeda's algorithm inputs, intermediate results and reconstructed wrapped phase surfaces. In, e. g. Figs. 5a and 5b there is the synthetically generated interferometric input image and its wrapped phase calculated directly from the disruption, in Figs. 5c and 5d – periodograms of interferogram and extracted carrier sideband, in Fig. 5e – the algorithm output: retrieved information about phase (in this case properly reconstructed) and in Fig. 5e – fringe pattern generated on the basis of the reconstructed wrapped phase surface.

In a typical case, low and intermediate (relatively to $f_0$) phase disruptions were reconstructed properly irrespective of fringe frequency (q. v. Fig. 5). For odd fringes count, there was a phase shift of $-\pi$, but further analysis showed that it was caused by the specificity of the fringe generator that was fixed to have zero phase and the local maximum of the cosine fringe function (i. e., $x=0$) in the image center. In the leftmost pixel, there was a half-period of the cosine function and the local minimum of the fringe value. When the FFT is evaluated, this local minimum is equivalent to $\pm\pi$ phase offset in the leftmost pixel associated with spatial point having $x=0$ coordinate.
<figure><img src="image_7.png" /><figcaption>Fig. 5. Takeda's algorithm with medium fringe frequency and moderate phase disruption</figcaption></figure>

Furthermore, for non-integer fringes count, there are some errors in the unwrapped phase that tend to accumulate near the border of the image, at the side from which the carrier side band was extracted (q. v. Fig. 6). The influence of above mentioned reconstruction faults is significant for small fringe frequencies and becomes almost negligible in the intermediate or higher bands. The significance of inter-pixel location of the true 1st order spectral peak is the most important cause of this effect, and the relative inaccuracy of peak localization, which in this case had 1px precision, is inversely proportional to the value of $f_0$. Some overcomes have been proposed to mitigate this undesired property, i. a. by determining Fourier peaks at sub-pixel level [13] (one attempt may be to infer about the peak exact location by utilizing the information of side band contour map). The applicability of these attempts is usually strongly restricted to interferograms with sufficiently low phase noise level and therefore the problem remains unsolved when the interferogram periodicity is strongly deteriorated (e. g. for tear film evaporating from contact lens surface).

Figures 7 and 8 depicts the conditions, in which proper phase reconstruction was impossible and the result was distorted by the excessive level of phase disruption. When there are closed loops in the interference pattern, it means that there is a region where decrease of instantaneous phase introduced by the disruption is too fast and cannot be compensated by the increase from constant frequency $f_0$ component. An area of local phase descent in one direction exhibits negative instantaneous frequency in that direction, as $f_x(x,y) = \frac{\partial \phi(x,y)}{\partial x}$. Negative frequency components from the $C(f-f_0,y)$ carrier
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Takeda's algorithm: non-integer fringe frequency, moderate phase disruption</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Takeda's algorithm: medium fringe frequency, intense phase disruption</figcaption></figure>

side band superimposes on the conjugate symmetric $C^*(f+f_0,y)$ side band negative frequency component and vice versa for positive frequencies. This introduces ambiguity in at least part of the area between two 1st order Fourier peaks. Thus, the side bands become partially indistinguishable and in general case it is impossible to eliminate the overlapping components from the chosen band. When the band is extracted using circular mask, a part of it still maintains information about its conjugate counterpart. IFFT of spectrum extracted in such manner amplifies or weakens some frequency components when compared to the unobtainable pure single band transformation. This eventually leads to asymmetry and some artefacts occurring in the reconstructed wrapped phase and also distorts the reconstructed fringe pattern. Although the visual disturbance is not necessarily of a substantial extent, in many cases it breaks some of the phase discontinuity loops or introduces open areas of such discontinuity. Such structures are associated with phase residues and the presence of such residues is a sign of improper phase surface reconstruction. What is more, it is impossible to remove phase discontinuities by unwrapping the phase along any close loop encompassing a single residue. This makes the phase unwrapping task non-trivial for residue-afflicted surfaces and its result will always be imperfect. Various phase unwrapping algorithms, however, can diminish residues influence to a different degree. Nonetheless, as a rule of thumb one can state that the more effective in accomplishing this task a method is, the more complicated its implementation would be and its computational complexity may be also disproportionately higher than the achieved improvement [4, 12].

In preceding examples involving synthetic fringe pattern, windowing was not applied, both in the spatial division (before FFT) and in the frequency division
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Takeda's algorithm: low (top row) and high (bottom row) fringe frequencies, intense phase disruption</figcaption></figure>

(before carrier side band masking) due to its negligible role as the perfect fringes periodicity and lack of spectral background noise were guaranteed.

**The Goldstein's Phase Unwrapping Algorithm.** Whereas phase unwrapping is trivial in 1-D problems, this is not true for higher dimensionality. In an ideal situation of perfectly reconstructed wrapped phase this would still be a relatively simple task. However, it is not the case applicable to most *in vivo* measurements. The phase in a complex image is subject to many factors that can disrupt its values. The influence of these phenomena reduces the useful signal detectability. Procedures have been developed that try to achieve a trade off between solution accuracy under these adverse conditions and computational requirements. No *standard procedure* to solve the phase unwrapping problem has been established, thus a large variety of algorithms are in use. Given the failure of a straightforward integration when applied to noisy data, other methods utilizing residue detection and branch cuts, quality maps or error minimization are commonly used. They can be divided into two main types [14]: a) path-following methods and b) minimum-norm methods. They differ in efficacy and computational complexity. Detailed explanation of their quality, robustness and execution time is presented in [4].

Goldstein's method (sometimes referred to as *Goldstein-Zebken-Werner algorithm*) has a fairly complex procedure structure [15]. Firstly, the polarity is computed for the phase residues. This involves phase gradient calculation, residues detection and residue charge determination (±1). Subsequently, the branch areas to cut are located, along which a phase integration path may not intersect. These are, i. a. lines linking residues of opposite polarity and, in some implementations, isolated residues located near image borders together with their connections. Typical LSI contact lens image produces a numerous amount of branch cuts as compared to the corneal inter-
ferograms, what indicates much higher level of image disruption caused by the drying
of the tear film. Finally, the instantaneous phase surface is reconstructed during flood
fill procedure that prevents every integration path from crossing the branch cuts.

*Advantages and Limitations of Goldstein's Algorithm.* Goldstein's algorithm takes a relatively short time to execute that makes one of the fastest phase unwrapping procedures. It also tries to reduce defects significance relative to a trivial phase unwrapping approach as depicted by Fig. 9 that shows failures during the unwrapping of imperfectly reconstructed phase. It takes measures to minimize the number of discontinuities and optimizes branch cuts location to locally minimize their length (q. v. Fig. 9b). Nonetheless, branch cuts placement can isolate segments of the image and thereby prevent them from being unwrapped. This occurs especially in datasets containing many residues. Figure 9b also reveals that the algorithm is subject to numerical pitfalls. In the utilized implementation, residues placed close to image borders are connected with them by two branch cuts directed at the angle of 45° instead of one at 90°. Thus a separated triangular areas are formed where no information about unwrapped phase was retrieved.

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** examples of phase unwrapping results (without phase residues masking). Discontinuities of the surface are present where the algorithms were unable to wrap phase due to residue-induced ambiguities</figcaption></figure>

There are some more efficient alternatives, but they are also of much higher
computational complexity [4, 16]. Much better (i. e. the most discontinuity-free)
solution is produced by, e. g. the minimum *l*^p*-norm algorithm. Average background
noise is also minimized in comparison to other algorithms. However, this is the second
slowest algorithm after the mask cut approach with execution time usually longer
than Goldstein in about one to two orders of magnitude (actual value vary greatly,
as it depends on the convergence speed [4]).

Some additional steps can be performed to locally improve the quality of this
method results. The modification of the Goldstein's algorithm described below al-
lows complete reduction of significant branch cuts at the expense of losing phase
information in some image parts outside its central area.

*Algorithm Modification.* After determining the residue polarities, the image is pre-
liminary masked to achieve a region, where the mask contains only regions with no
residues. Also accepted are the areas with pairs opposite charge residues separated
by a short distance. All subsequent actions may be performed only within the masked
region. This condition ensures lack of any significant discontinuities in unwrapped
phase. The mask is calculated as a circle expanding (by default from the image center)
as long as it fulfils the above conditions. Then, four more circles are created on its
border and expanded. This subroutine is performed in a recursive manner and thus
the modification is called a *recursive bubble* approach.

## 2.2 Phase Parameters as TFSQ Descriptors

Three phase components were derived from the phase coordinates (two spatial and one
intensity component). The variance of all of them was significantly different than
the variance of the phase itself. However, after applying the linear detrending of the un-
wrapped phase, the variance of phase was much smaller than of spatial components.
Thus, PCA combined and transformed the coordinates into first two principal com-
ponents, whilst the intensity of the phase map was assigned virtually intact to last
component. Therefore phase and 3rd component variations were the same. This
indicates that PCA procedure is redundant and phase variance is equally adequate
TFSQ descriptor.

# 3 Results

There was an application developed that underwent tests on a set of static and dynamic
data. The initial input data set consisted of 19 static PNG images and 13 AVI video se-
quences. Results for the synthetic undistorted images were consistent with theoretically
predicted results. The output for *in vivo* LSI (q. v. Fig. 10) has shown that the al-
gorithm is capable of reconstructing phase for slightly distorted images, but gives
ambiguous output in case of strongly deteriorated film. Processing of a single movie in
MATLAB environment lasted between 30 min and 16 h depending on disorder level.

4 Discussion and Conclusions

Spectral representation of LSI images provide mostly information of sub-global char-
acter related to the frequency of fringes. On the other hand, the phase extracted
from those images carries information of local character. The unwrapped phase has
no simple physical interpretation. Although it is not directly related to the wavefront,
it is causally associated with it. However, there are similarities between the two con-
sidered representations. In case of deteriorating quality of tear film covering the eye
cornea, the spectral representation of LSI images becomes less regular and the phase
representation starts to lose its continuity.
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** exemplary outputs of the tear film surface quality assessing system</figcaption></figure>

**Conclusions** Analysis of the measurements and theoretical background suggests that the method joining Takeda's and Goldstein's algorithms with phase variability analysis as tear film surface quality descriptor is capable of supplying useful information. Still, along with the increase of the phase noise occurring when tear film dries up, the efficacy of this solution decreases significantly.

Tests on LSI images, frequently exhibited saddle-shaped formation in phase surface and chessboard-like structures in the phase maps. They were consistent with the visible input phase variations. Their occurrence might arise due to optical aberrations of objective lens and the natural corneal asphericity. Thus the method has a potential for serving in other closely-related applications, like testing and calibration of medical diagnostic equipment.

*Future Recommendations.* As a way of overcoming issues encountered, a enhanced or new method of interferogram analysis may be adapted. In particular an extended set of both wrapped and unwrapped phase parameters might be considered, e. g.: *a*) excess kurtosis; *b*) texture statistics derived from gray-level co-occurrence matrix (GLCM); homogeneity, energy/uniformity, correlation and contrast; *c*) Rényi entropy and *d*) Kass-Witkin coherence [17]. Thereafter a modification of hardware in experimental setup can be taken into account to improve the input data quality.

## Acknowledgments

The authors would like to express their gratitude to **Prof. Ewa Skubalska-Rafajłowicz** and **Dr. Dorota Szczęsna-Iskander** for extensive and fruitful cooperation and for their encouragement, advice and support.

Supported by the European Regional Development Fund within Innovative Economy Operational Program co-financed by the Foundation for Polish Science (POMOST/2012-5/8/0072).
5 References

[1] D. H. Szczęsna, D. Alonso-Caneiro, D. R. Iskander, S. A. Read, and M. J. Collins, "Lateral shearing interferometry, dynamic wavefront sensing, and high-speed videokera-toscopy for noninvasive assessment of tear film surface characteristics: A comparative study", *J. Biomed. Opt.*, vol. 15, no. 3, 2009, pp. 1-9.

[2] D. H. Szczęsna and D. R. Iskander, "Robust estimation of tear film surface quality in lateral shearing interferometry", *J. Biomed. Opt.*, vol. 14, no. 6, pp. 064039-1-064039-4, 2009.

[3] M. Takeda, H. Ina, and S. Kobayashi, "Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry", *J. Opt. Soc. Am.*, vol. 72, no. 1, pp. 156-160, Jan. 1982.

[4] S. Braganza, *Analysis of phase-unwrap quality of OQM images*, URL: http://www.coe.neu.edu/Research/rcl/projects/phaseunwrap/Analysisofphase-ver2.htm, May 2006.

[5] M. Rottenkolber and H. Podbielska, "High precision Twyman-Green interferometer for the measurement of ophthalmic surfaces", *Acta Ophthalmol. Scand.*, vol. 74, no. 4, pp. 348-353, 1996.

[6] T. J. Licznerski, H. T. Kasprzak, and W. Kowalik, "Application of Twyman-Green interferometer for evaluation of in vivo breakup characteristic of the human tear film", *J. Biomed. Opt.*, vol. 4, no. 1, pp. 176-182, 1999.

[7] D. Malacara, M. Servín, and Z. Malacara, *Interferogram analysis for optical testing*, 2nd ed., Taylor & Francis Group, Ed., ser. Optical engineering (Marcel Dekker, Inc.) Boca Raton, FL: CRC press, 2005, vol. 84.

[8] D. H. Szczęsna, J. Jaroński, H. T. Kasprzak, and U. Stenevi, "Interferometric measure-ments of dynamic changes of tear film", *J. Biomed. Opt.*, vol. 11, no. 3, pp. 034028-1-034028-8, 2006.

[9] T. J. Licznerski, "Interferencyjne i modelowe badanie filmu łzewowego oka ludzkiego", Polish, PhD thesis, Politechnika Wrocławska, Wrocław, 1998.

[10] M. Takeda and K. Muto, "Fourier transform profilometry for the automatic mea-surement of 3-D object shapes", *Appl. Opt.*, vol. 22, no. 24, pp. 3977-3982, Dec. 1983.

[11] R. M. Goldstein, H. A. Zebker, and C. L. Werner, "Satellite radar interferometry: Two-dimensional phase unwrapping", *Radio Sci.*, vol. 23, no. 4, pp. 713-720, 1988.

[12] D. C. Ghiglia and M. D. Pritt, *Two-dimensional phase unwrapping: Theory, algorithms, and software*, ser. Wiley-Interscience. Wiley, 1998.

[13] E. Mainsah and D. G. Greenwood Chetwynd, *Metrology and Properties of Engineering Surfaces*, E. Mainsah, J. A. Greenwood, and D. G. Chetwynd, Eds., ser. Optoelectronics, Imaging and Sensing. Springer, 2001.

[14] R. Gens, *Phase unwrapping*, URL: http://www.asf.alaska.edu/~rgens/teaching/geos639/phase_unwrapping.pdf, 2006.

[15] R. M. Goldstein and C. L. Werner, "Radar ice motion interferometry", in *Proceedings of 3rd European Remote Sensing Symposium "Space at the Service of our Environment"*, vol. 2, Florence, Italy, 1997, pp. 969-972.

[16] R. Gens, "Two-dimensional phase unwrapping for radar interferometry, Developments and new challenges", *International Journal of Remote Sensing*, vol. 24, no. 4, pp. 703-710, 2003.

[17] M. Kass and A. Witkin, "Analyzing oriented patterns", *Computer vision, graphics, and image processing*, vol. 37, no. 3, pp. 362-385, 1987.