Method of surface topography retrieval by direct solution of sparse weighted seminormal equations

Jeffrey Koskulics; Steven Englehardt; Steven Long; Yongxiang Hu; Knut Stamnes

doi:10.1364/OE.20.001714

Introduction

Surface topography is often characterized indirectly by measurements of slopes. An interesting class of examples can be drawn from instruments that make high resolution measurements of wind roughened water surfaces [1–6], whose output data is in the form of slopes. Such measurements are important in the continuing investigation of rough ocean surfaces [7], whose characteristics play a central role in models of radar scatterometry [8], lidar [9] and sunglint [10]. It has been shown that accurate models of radiative transfer at the ocean surface for near infrared and shorter wavelength bands require characterization of topographic details beyond simple slope distributions [11]. Topography data can yield significant new insights into the understanding of the rough ocean surface, but enhanced understanding requires some method of estimating topography from slope data as no other direct method of detailed topographic measurement exists.

Examples of topographic recovery techniques include approaches based on integration [12–14], and expansion of solutions using orthogonal functions, such as sinusoids (Fourier transforms) [15, 16], radial basis functions [17] and Legendre polynomials [18] using slope data as input. Topography can also be recovered using measurements of curvature [18, 19]. These methods can suffer from drawbacks including artifacts caused by the global propagation of errors for path integration approaches, sensitivity to high frequency noise for Fourier transforms, or poor numerical performance due to the high order polynomials that are needed to treat highly curved rough surfaces. Furthermore, it is not clear how a priori data, such as the mean and variance of elevation, slope and curvature, can be considered alongside measurements in these methods.

The problem of surface topography retrieval has many similarities to phase retrieval in adaptive optics (AO) systems. In many AO systems, the slopes of a distorted incoming plane wave are measured, and a compensation is determined and applied, such as by deforming a mirror into a shape that matches the distorted wave front. The term “wave front phase” in the cited literature is analogous to the surface topography considered here. An important class of phase retrieval methods is based on least squares solution of a coupled system of finite difference equations relating wave front phase and slope [20–22]. The equations can be formulated as a matrix-vector multiplication [23] whose form is similar to the sparse system of equations encountered in the theory of networks [24], which can be solved by sparse equation solvers. In this way, large systems can be stored in memory of practical size and solved using iterative techniques, such as the method of successive over relaxation (SOR) [25].

In a variety of physical situations, a sparse system of coupled linear equations results from finite difference equations, such as those used to represent measured gradient fields [26], and this approach can be generalized to include other quantities such as surface elevation and curvature. A priori constraints are often a practical necessity to ensure full system rank. Immediate measurements and a priori observations having different units and precision can be weighted so that each quantity contributes appropriate constraints to the solution.

Until recently, computer memory restrictions and computational costs have inhibited attempts at direct solution of very large systems of equations containing millions of state and measurement vector elements, such as those resulting from megapixel slope image data. Advances in computational hardware and software, including the introduction of 64-bit operating systems, multithreaded/multicore microprocessor architectures, large low-cost memory arrays and powerful new solution algorithms make possible direct solution of large systems of coupled linear algebraic equations.

Developed in this paper is a method of surface topography retrieval via sparse weighted linear least squares. The retrieval problem is formulated in such a way that the solution of the inverse problem allows for simultaneous consideration of direct measurements of surface slopes, a priori observations such as mean surface curvature and elevation without it being necessary to directly compute the inverse of the system. This new formulation and the corresponding method of solution can be applied to solve a wide variety of classical measurement problems, examples of which include topographic retrieval from two dimensional arrays of one-dimensional slope measurements [1, 2], two-dimensional arrays of two-dimensional slope data [3, 4, 27], and where individual measurement data are irregularly sampled, missing or have varying uncertainty [1]. As such, it can be viewed as a generic methodology that can be adjusted to suit the particular retrieval problem of interest. To demonstrate the power and practical utility of the method, examples are provided in which it is applied to large linear systems based on synthetic data as well as real data.

Mathematical formulation

Surface elevation and slopes

A surface can be approximately described by a discrete set of vertices (points on the surface) arranged as a mesh to form a collection of triangular planar facets. The slope of a planar facet can be characterized by a surface normal vector, computed by the vector cross product shown in Fig. 1 .

Fig. 1 A triangular planar facet (hatched area) is bounded by corners having position vectors ${\vec{r}}_{a}$ , ${\vec{r}}_{b}$ and ${\vec{r}}_{c}$ in a Cartesian coordinate system. The surface normal vector $\vec{n}$ can be computed by cross product and normalized to unit length by $\hat{n} = \vec{n} / | \vec{n} |$ .

Download Full Size | PDF

Using the notation ${\vec{r}}_{a} = [\begin{matrix} x_{a} & y_{a} & z_{a} \end{matrix}]$ for the Cartesian $x$ , $y$ and $z$ components of the vertex position vector, the surface normal can be computed by the cross product $\vec{n} = ({\vec{r}}_{b} - {\vec{r}}_{a}) \times ({\vec{r}}_{c} - {\vec{r}}_{a})$ , which expands to:

\vec{n} = \begin{matrix} [(y_{b} - y_{a}) (z_{c} - z_{a}) - (z_{b} - z_{a}) (y_{c} - y_{a})] \hat{i} + \\ [(z_{b} - z_{a}) (x_{c} - x_{a}) - (x_{b} - x_{a}) (z_{c} - z_{a})] \hat{j} + \\ [(x_{b} - x_{a}) (y_{c} - y_{a}) - (y_{b} - y_{a}) (x_{c} - x_{a})] \hat{k} \end{matrix} .

where

\hat{i}

,

\hat{j}

and

\hat{k}

are Cartesian unit vectors pointing along the

x

,

y

and

z

axes of the coordinate system, respectively.

The horizontal ( $x$ and $y$ ) components of the vertex positions, and their differences, can be treated as the constant parameters: $c_{1} = x_{b} - x_{a}$ , $c_{2} = x_{c} - x_{a}$ , $c_{3} = y_{b} - y_{a}$ and $c_{4} = y_{c} - y_{a}$ . The cross product simplifies to:

\begin{matrix} [c_{3} (z_{c} - z_{a}) - c_{4} (z_{b} - z_{a})] \hat{i} + \\ \vec{n} = [c_{2} (z_{b} - z_{a}) - c_{1} (z_{c} - z_{a})] \hat{j} + \\ [c_{1} c_{4} - c_{2} c_{3}] \hat{k} . \end{matrix}

Surface normal unit vectors can be formed by scaling the surface normal vector of Eq. (2) through division by its magnitude $\hat{n} = \vec{n} / | \vec{n} |$ , where $| \vec{n} | = {(\vec{n} \cdot \vec{n})}^{1 / 2}$ . Because the vertical component of the surface normal vector is equal to the constant $c_{1} c_{4} - c_{2} c_{3}$ , its magnitude can be recovered from a surface normal unit vector by equating vertical components, $\vec{n} \cdot \hat{k} = | \vec{n} | \hat{n} \cdot \hat{k}$ . Surface normal vectors given with unit length can be rescaled to give Eq. (2) by multiplication with the magnitude of the surface normal vector, which equals:

| \vec{n} | = (c_{1} c_{4} - c_{2} c_{3}) / \hat{n} \cdot \hat{k} .

Arbitrary arrangements of surface normal data provided on non-uniform horizontal grids can be related to vertex elevations using Eqs. (2, 3). In many cases, slopes are provided on a grid that has a uniform spacing so that horizontal vertex coordinates can be organized in perpendicular rows and columns. This arrangement enables further simplification by aligning the $x$ and $y$ coordinate axes with the rows and columns of vertices, so that $x_{a} = x_{c}$ and $y_{a} = y_{b}$ . In this case, $c_{2} = c_{3} = 0$ . Surface normal vector components become functions of just two vertices:

\vec{n} \cdot \hat{i} = c_{4} (z_{a} - z_{b})

and

\vec{n} \cdot \hat{j} = c_{1} (z_{a} - z_{c}) .

For simplicity, the special case of row-column alignment with the $x$ and $y$ coordinate axes will be considered throughout the remainder of this paper, though this method can be applied to arbitrary grids of slopes without loss of generality.

Arrays of slope data

In this formalism, an approximation of the surface is represented by a discrete collection of vertices, and gradient measurements are represented by surface normal vectors that are perpendicular to triangular planar surface facets. In slope imaging systems, pixel data represent measurements of the local surface gradient. Each pixel has a finite area with a well-defined outline in the object plane that suggests a choice of mesh geometry placing the horizontal position of vertices at the pixel corners.

The vertices and surface normals can be ordered by column, an example of which is illustrated in Fig. 2 . Because a plane is well-defined by three distinct points, slope measurements are considered to occur on triangular planar facets. In this configuration, one planar facet corresponds to each pixel, and the triangular facet in the lower right corner of each pixel is effectively ignored.

Fig. 2 Example mesh shows the ${\hat{n}}_{s}$ surface normal unit vector measurements, and the $z_{t}$ vertex elevations. The shaded areas represent locally planar facets on which the measurements are considered to occur.

Download Full Size | PDF

In this example, there are $M = 4$ rows and $N = 3$ columns of surface normal measurements, and $M + 1 = 5$ rows and $N + 1 = 4$ columns of vertices. Surface normal vectors and vertices residing at row $m$ and column $n$ are numbered with subscripts $s = m + M (n - 1)$ and $t = m + (M + 1) (n - 1)$ , respectively.

Matrix-vector formulation

The surface normal vector measurements can be grouped by Cartesian component to give two measurement vectors:

b_{x} = {[\begin{matrix} {\vec{n}}_{1} \cdot \hat{i} & {\vec{n}}_{2} \cdot \hat{i} & \dots & {\vec{n}}_{s} \cdot \hat{i} & \dots & {\vec{n}}_{M \times N} \cdot \hat{i} \end{matrix}]}^{T}

and

b_{y} = {[\begin{matrix} {\vec{n}}_{1} \cdot \hat{j} & {\vec{n}}_{2} \cdot \hat{j} & \dots & {\vec{n}}_{s} \cdot \hat{j} & \dots & {\vec{n}}_{M \times N} \cdot \hat{j} \end{matrix}]}^{T}

for the horizontal and vertical components, respectively (the superscript

^{T}

indicates vector transpose). Vertex elevations are similarly grouped in a column state vector:

z = {[\begin{matrix} z_{1} & z_{2} & \dots & z_{t} & \dots & z_{(M + 1) (N + 1) - 1} \end{matrix}]}^{T}

Related by finite difference equations, the measurement vectors $b_{x}$ and $b_{y}$ are linear functions of the state vector $z$ and can be expressed as matrix-vector products: $b_{x} = A_{x} z$ and $b_{y} = A_{y} z$ , constructed from Eqs. (4, 5), respectively. Row and column numbering is determined by the mesh numbering system, such as in the example of Fig. 2.

Rank

A system of equations containing only slope expressions is rank deficient due to an arbitrary constant, the so-called “piston” term, that can be added to a solution without changing slopes [24]. To resolve this ambiguity, a single equation of average elevation can be added to the system to restore full rank.

Rank deficiency can also arise from “island” vertices or groups of vertices that share no adjacency. Such deficiency exists in the example of Fig. 2 where the 20^th vertex is unconnected by surface normal data to the rest of the mesh. Islands can also arise where one or more groups of missing or unreliable data isolate groups from each other. A single expression of average elevation generally does not solve the problem of rank deficiency arising from topographic islands because an arbitrary constant can be added to one island and subtracted from another without affecting the average elevation.

A priori constraints

To ensure full system rank, the inclusion of additional relationships is necessary and can be accomplished in a variety of ways, such as by direct measurements of the elevation of one or more vertices. Alternatively, a priori observations may be added. For example, if the mean elevation is known, each vertex can be constrained by augmenting the system with an identity matrix that multiplies the state vector and set equal to the mean elevation: $I z = \bar{z}$ , where each row added has the form $z_{t} = \bar{z}$ .

Smooth solutions can be constrained on the basis of a priori curvature. Curvature can also be considered as immediately measured data, such as in the case of curvature measurements [19]. The curvature of plane curves that are specified explicitly as a function of one horizontal coordinate can be approximated by a second derivative $κ = z^{″} {(1 + {z^{'}}^{2})}^{- 3 / 2} \approx z^{″}$ for $z^{'} \approx 0$ . By choosing planes parallel to the rows and columns (where one horizontal coordinate is held constant), the second derivative can be approximated by the finite difference equations:

κ_{x}_{t} \approx c_{1}^{- 2} (z_{m - 1, n} - 2 z_{m, n} + z_{m + 1, n})

and

κ_{y}_{t} \approx c_{4}^{- 2} (z_{m, n - 1} - 2 z_{m, n} + z_{m, n + 1}),

where the subscripted x and y indicate curvature along columns and rows, respectively. Systems of equations

K_{x} z = κ_{x}

and

Κ_{y} z = κ_{y}

can be constructed from Eqs. (9) and (10), respectively, and set equal to the mean curvature (usually zero):

{\bar{κ}}_{x} = {\bar{κ}}_{y} = 0

. A system of measured x- and y-slopes, with a priori x- and y-plane curvatures and elevation can be combined in a block matrix

A

:

A z = [\begin{matrix} A_{x} \\ A_{y} \\ Κ_{x} \\ Κ_{y} \\ I \end{matrix}] z = [\begin{matrix} b_{x} \\ b_{y} \\ κ_{x} = 0 \\ κ_{y} = 0 \\ z = \bar{z} \end{matrix}]

In general, there is considerable flexibility in choosing which measurements and a priori constraints to include. Blocks can be added or deleted as needed depending on the problem at hand. For example, arrays of one-dimensional slope data [1, 2] could be treated by setting $b_{x}$ equal to the slope measurements. If average slope in the perpendicular direction were known to be zero, $b_{y}$ could be set to zero, or $A_{y}$ and $b_{y}$ could be deleted from the system altogether.

An important restriction worth mentioning is that for large systems, memory considerations require that any added equation be sparse, generally limiting systems to those based on finite difference equations or identity relations. The structure of a sparse system can be inspected by plotting the non-zero row-column entries. Figure 3 shows the non zero entries of the sparse matrix given by the linear system of Eq. (11) using the example shown in Fig. 2. In the following sections, Eq. (11) will be considered for further analysis using synthetic and real data.

Fig. 3 The non zero entries of the sparse matrix given by the linear system of Eq. (11) using the example of Fig. 2.

Download Full Size | PDF

Weighted least squares

_{$A z = b$} has a greater number of rows than columns which results in an over-determined system. Over-determined systems generally have no exact solution $z$ . Instead, an approximate least squares solution is computed to minimize the 2-norm of a residual error vector ${‖ A z - b ‖}^{2} = {(A z - b)}^{T} (A z - b)$ by solving the normal equations: $A^{T} A z = A^{T} b$ .

Biases caused by differences in the relative uncertainty of measurement components can be compensated by weighting the individual components of the residual error vector by a weighting matrix $W$ giving: ${‖ W (A z - b) ‖}^{2} = {(A z - b)}^{T} C (A z - b)$ , where $C = W^{T} W$ . The resulting equations minimize the weighted error:

A^{T} C A z = A^{T} C b .

_$C$ can be specified by a covariance matrix $Σ$ , where each matrix element specifies an estimate of the error covariance between measurement components: $Σ_{i j} = cov (b_{i}, b_{j})$ . The choice of $C = Σ^{- 1}$ leads to the Best Unbiased Linear Estimate (BLUE) as stated by the Gauss-Markov theorem [28] which gives the lowest mean squared error of the estimate. Specification of variance can be accomplished in several ways: variance of immediately measured components can be determined from repeated measurements of an unchanging state, and a priori variance is frequently specified explicitly as a parameter of a statistical distribution.

It should be mentioned that for surfaces that are relatively smooth on small scales, the elevation of neighboring vertices is likely to have some non-zero covariance that decreases with increasing horizontal separation. Therefore, off-diagonal terms will appear in $Σ$ . Generally, inversion of an off-diagonal $Σ$ will fill in nonzero elements so that the resulting matrix $C$ is full and often much larger than the available memory. This practical limitation generally prevents the consideration of off-diagonal covariance in large systems, despite any additional information that a full covariance matrix might contain. In the following analysis, the variance matrix $Σ_{i i} = var (b_{i})$ is used instead and off-diagonal terms are ignored ( $Σ_{i \neq j} = 0$ ).

Solution method

Sparse orthogonal-triangular (QR) factorization is an efficient, numerically stable method of solution for overdetermined systems [29], supported in the backslash “\” operator and ^qr and ^lscov functions in recent versions of Matlab subsequent to R2009b, also available as C + + source code in the SuiteSparse archive [30]. QR factorization expresses a matrix as the product of two matrices $A = Q R$ , where $R$ is upper-triangular, and $Q$ is an orthogonal matrix, such that $Q^{- 1} = Q^{T}$ and $Q^{T} Q = I$ . In general, $R$ is sparse and $Q$ tends to be full.

Column reordering is essential to maximize the sparsity of the $R$ factor (for efficient storage in memory) and to maintain practical limits on factorization time. Minimum degree column reordering algorithms aim to reduce the number of nonzero matrix elements by maximizing the number of elimination steps that can be skipped in the factorization process. Column reordering replaces $A$ and $z$ in Eq. (12) with the column-permuted matrix $A_{p} = A P$ and the row-permuted vector $z_{p} = P^{- 1} z$ , where $P$ is an orthogonal permutation matrix, where $P^{- 1} = P^{T}$ and $P^{T} P = I$ .

The normal equations can be solved by factoring the matrix $W A_{p} = Q R$ . In cases where storage of the factor $Q$ would exceed available memory resources, the $Q$ factor is dropped, and the factorization is referred to as “Q-less” QR factorization. Substituting these factors into Eq. (12) gives the seminormal equations [31]: ${(Q R)}^{T} Q R z_{p} = {(A_{p})}^{T} C b$ . The left side simplifies to $R^{T} R z_{p}$ , and the right side of the seminormal equation can be rewritten as the vector $u = {(A_{p})}^{T} C b$ , giving:

R^{T} R z_{p} = u .

The advantage of solving the seminormal equations by factorization is that $A$ can be reordered and factored once and $R$ can be repeatedly applied to the solution of multiple measurement vectors, not all of which may be known or can be stored in memory at the time of factorization. The product of two triangular matrices is efficiently solved by forward and back substitution for $z_{p}$ . Errors that accumulate due to limited numerical precision can be reduced by solving the corrected seminormal equations [32]. An error term is computed by substituting the solution $z_{p}$ into Eq. (13), $r = W (A_{p} z_{p} - b)$ . This error can be used to compute a correction factor by solving $R^{T} R Δ z_{p} = {(A_{p})}^{T} C r$ for a correction factor $Δ z$ that can be added to find a new solution of reduced error: ${z^{'}}_{p} = z_{p} + Δ z_{p}$ . After the corrected solution is computed, a least squares solution with proper ordering can be restored by multiplication with the permutation matrix: $z = P {z^{'}}_{p}$ .

The Matlab backslash operator can be used to perform all of these steps and solve the weighted normal equations $W A z = W b$ using the syntax: z = (W*A)\(W*b). Multiple measurement vectors $b_{1}$ , $b_{2}$ … $b_{n}$ can be combined for efficient solution by horizontal concatenation into a matrix $B$ , where the columns of $B = [\begin{matrix} b_{1} & b_{2} & \dots & b_{n} \end{matrix}]$ . In this way, Matlab computes the intermediate factorization step only once and applies the results to all of the measurement vectors. Solution vectors appear as columns of a solution matrix $Z$ .

Synthetic data retrieval

The method of topographic reconstruction is illustrated in an example using synthetic data. While no single example will capture the full range of possible topographies that may be encountered in real data, much can be learned by examining synthetic data. The procedure will be as follows:

1. A surface is given by an elevation function of two horizontal coordinates
2. A set of synthetic slope measurements is created by evaluating analytical derivatives of the elevation function on a grid of points and adding random numbers to the results to simulate noise
3. Weights for slope and a priori measurement components are chosen
4. Surface elevation is estimated by the procedure described above.

The surface chosen is the analytical function:

z (x, y) = \sin (x / 3) \sin (y)

with analytical derivatives:

\frac{\partial z}{\partial x} = \cos (x / 3) \sin (y) / 3

and

\frac{\partial z}{\partial y} = \sin (x / 3) \cos (y) .

An ideal surface is constructed by evaluating Eq. (14) on a grid of 641x481 points over a spatial domain that shows a handful of peaks. Synthetic slope measurements are computed on a grid of 640x480 points by evaluating Eq. (15) halfway between the columns, and by evaluating Eq. (16) halfway between the rows of the grid used to evaluate the ideal surface.

A priori variance is computed by evaluating Eqs. (14–16) and the second derivatives of Eq. (14) (for curvature) on a grid of points, followed by computing the variance of the resulting sets. Two cases using synthetic slope measurements are made by contaminating the synthetic slope measurements with artificial noise by adding normally distributed random numbers with specified standard deviation. The weighting of slope measurements is determined by the standard deviation of the added noise.

A low noise set of synthetic x- and y-slope measurements data is first examined. Figure 4 shows a retrieval based on nearly ideal synthetic data, where the added noise has a standard deviation of 5.118 $\times$ 10⁻⁴, one one-thousandth of the a priori y-slope data standard deviation (a zero-noise case would have infinite weighting coefficients and cannot be analyzed by this method). The retrieved surface appears to be smooth. The difference between the original and retrieved surfaces shows deviations on the order of 10⁻⁴ times the maximum elevation with a standard deviation of error 5.2687 $\times$ 10⁻⁵. The errors appear largely random, though signs of a periodic artifact reminiscent of the original synthetic data can be seen in this difference.

Fig. 4 Left, surface retrieved from low-noise slope data. Right, retrieval error as the difference between retrieved and ideal surface. Surface color is a function of elevation shown in the color bars.

Download Full Size | PDF

Noisy slope data with noise having a standard deviation of 0.0512 (which is 30% and 10% of the x- and y-slope a priori standard deviation, respectively) results in retrievals with more distortion, shown in Fig. 5 . The distortion appears smooth, with the largest errors less than 2% of the original surface’s maximum elevation appearing at the boundaries of the solution. Perhaps the errors in slope balance each other in the interior and errors caused by the elevation constraints become apparent at the boundaries. Most importantly, the standard deviation of error in elevation is 0.0055, or about 0.5% of the original surface’s maximum elevation.

Fig. 5 Left, surface retrieved from noisy slope data. Right, solution error as the difference between retrieved and ideal surface.

Download Full Size | PDF

These computations were performed with Matlab (R2010b) using an Intel Core i3 processor on a personal laptop computer with 8 GB memory. Column reordering and factorization required 25 seconds with the resulting factor residing in 632 MB stored as a double precision sparse matrix. Solution of the column reordered prefactored corrected seminormal equations computed by forward and backward substitution using the backslash operator takes 2.86 seconds of computational time.

Water topography retrieval

Slope data were obtained from a refractive color slope imaging system with resolution 1024x1280 pixels in a construction similar to the surface gradient detector described by Zhang and Cox [3]. Imaging data were collected during an experiment in a wind wave tank at the NASA Air-Sea Interaction Research Facility at Wallops Island, Virginia with an example shown in Fig. 6 . The tank geometry is approximately 18x1x1.2 meter in length, width and height, respectively. The water fills approximately 2/3 of the vertical space. Wind blows along the longest dimension of the tank. Image data were collected approximately 9 meters from the air inlet. Hence, the waves are fetch-limited, just beginning the growth process of increasing in amplitude and wavelength.

Fig. 6 Left: color image data from the refractive slope imaging system. Middle and right, x and y surface normal unit vector components, respectively.

Download Full Size | PDF

The measurement system also employed capacitive elevation sensing probes that are visible in the corners of the image as dark pixels. These pixel data are considered to be contaminated so that no slope information can be obtained. Contaminated pixels are identified and set to zero slope before further analysis. Horizontal components of the surface normal unit vector were retrieved from the color image data using a newly developed method that will be detailed in an upcoming publication.

The variance matrix was scaled based on the height, slopes, and curvature having standard deviations of 0.1 meter, 0.058 meters per meter, and 500 per meter, respectively. The elevation standard deviation was set much greater than the observed wave amplitude, acting only loosely as a constraint. The standard deviation of slope was estimated directly from repeated measurements of a flat surface, corresponding to a one-sigma slope error of approximately 3.3° degrees, which may be considered relatively noisy slope measurements. Curvature standard deviation was established by comparing repeated solutions having different curvature standard deviations. High values of standard deviation had no appreciable effect on smoothness, whereas very low values tended to “blur” the surface, obscuring the smallest waves entirely. A compromise was established using the standard deviation at which a small smoothing effect became noticeable.

The retrieval results shown in Fig. 7 demonstrate well-behaved topographical solutions. All topographic points are found very close to the mean elevation, qualitatively similar to how wind waves appear by eye. The amplitudes of the retrieved wave crests are approximately the same amplitude as observed through the glass walls of the tank. Areas contaminated by identifiably bad data do introduce observable artifacts. In these areas, isolated islands of vertices that are unconnected by slopes can be found. Here, a priori data take over to fill in the blanks with reasonable estimates. Overall, the retrieved topography is free of obvious defects, which indicates that the solution is reasonable. Future work will examine in detail errors in slope and elevation using a standard reference surface with a shape that can be measured independently.

Fig. 7 Single frame example of retrieved elevation rendered using an artificial illumination model (Media 1).

Download Full Size | PDF

In terms of computational performance, the factorization procedure required 3.4 GB of memory and 95 seconds of computational time using an Intel Core i7 processor on a personal desktop computer. Prefactored seminormal equations are computed by elimination using 5.2 seconds of computation time and requires an additional 3.4 GB of memory as two copies of the factor are held in memory during the solution process.

Conclusions

A novel approach to the retrieval of surface topography is formulated in which the solution of the inverse problem allows for simultaneous consideration of immediately measured and a priori observations. Observations that can be included in this framework are modeled as coupled finite difference equations that result in a sparse linear system of algebraic equations. Observables such as elevation, slope and curvature can be included, treated as either immediate measurements or a priori constraints. By including such a priori observations, full system rank can be ensured, despite noisy or missing data points. Solutions are computed directly by multifrontal QR factorization without inversion of the full system matrix.

Synthetic data as well as imaging slope data obtained in a wave tank are used to demonstrate that solutions can be computed using least squares by combining noisy slope measurements with a priori elevation and curvature data, with each observation weighted by estimates of variance. Solutions retrieved from synthetic data have error standard deviations of approximately 1/10th of slope measurement error standard deviation. Solutions retrieved from imaging data are well behaved even in places where pixel data are known to be corrupt. The computations for the examples shown were achieved using a personal laptop and desktop computer. Reasonable solution times enabled the analysis of long time series of wind wave slope images, an example of which is shown in the companion movie file.

The enabling factors for this new method are computing resources of modest cost that have recently become available. The equations can be solved efficiently by utilizing capabilities of multicore/multithread processing and large memory arrays. The examples used to illustrate the method can be considered as members of a broader class of least squares retrievals that are also treatable by this method. This new formulation and the corresponding method of solution can be applied to solve a wide variety of other classical measurement problems. It is expected that these new developments will open new avenues to develop solutions for topographic and other more general sparse retrieval problems.

Acknowledgments

The authors acknowledge support from the NASA Calipso Project and Radiation Science Program. The authors would also like to thank the anonymous reviewers for their detailed comments and suggestions. S. Long wishes to thank NASA HQ for their support over the past 36 years, right up to this final work that marks the facility’s end of operations.

References and links:

1. C. Cox and X. Zhang, “Contours of slopes of a rippled water surface,” Opt. Express 19(20), 18789–18794 (2011). [CrossRef] [PubMed]

2. W. C. Keller and B. L. Gotwols, “Two-dimensional optical measurement of wave slope,” Appl. Opt. 22(22), 3476–3478 (1983). [CrossRef] [PubMed]

3. X. Zhang and C. S. Cox, “Measuring the two-dimensional structure of a wavy water surface optically: a surface gradient detector,” Exp. Fluids 17(4), 225–237 (1994). [CrossRef]

4. C. J. Zappa, M. L. Banner, H. Schultz, A. Corrada-Emmanuel, L. B. Wolff, and J. Yalcin, “Retrieval of short ocean wave slope using polarimetric imaging,” Meas. Sci. Technol. 19(5), 05503 (2008). [CrossRef]

5. Q. Li, M. Zhao, S. Tang, S. Sun, and J. Wu, “Two-dimensional scanning laser slope gauge: measurements of ocean-ripple structures,” Appl. Opt. 32(24), 4590–4597 (1993). [CrossRef] [PubMed]

6. B. Jähne, J. Klinke, and S. Waas, “Imaging of short ocean wind waves: a critical theoretical review,” J. Opt. Soc. Am. A 11(8), 2197–2209 (1994). [CrossRef]

7. W. Munk, “An inconvenient sea truth: spread, steepness, and skewness of surface slopes,” Annu. Rev. Mar. Sci. 1(1), 377–415 (2009). [CrossRef] [PubMed]

8. A. Freedman, D. McWatters, and M. Spencer, “The Aquarius scatterometer: an active system for measuring surface roughness for sea-surface brightness temperature correction,” presented at the IEEE International Geoscience & Remote Sensing Symposium, Denver, Colorado, July 31- August 4, 2006.

9. Y. Hu, K. Stamnes, M. Vaughan, J. Pelon, C. Weimer, D. Wu, M. Cisewski, W. Sun, P. Yang, B. Lin, A. Omar, D. Flittner, C. Hostetler, C. Trepte, D. Winker, G. Gibson, and M. Santa-Maria, “Sea surface wind speed estimation from space-based lidar measurements,” Atmos. Chem. Phys. Discuss. 8(1), 2771–2793 (2008). [CrossRef]

10. F. M. Bréon and N. Henriot, “Spaceborne observations of ocean glint reflectance and modeling of wave slope distributions,” J. Geophys. Res. 111(C6), C06005 (2006). [CrossRef] [PubMed]

11. S. Kay, J. Hedley, S. Lavender, and A. Nimmo-Smith, “Light transfer at the ocean surface modeled using high resolution sea surface realizations,” Opt. Express 19(7), 6493–6504 (2011). [CrossRef] [PubMed]

12. R. Klette and K. Schlüns, “Height data from gradient fields,” Proc. SPIE 2908, 204–215 (1996). [CrossRef]

13. K. Schlüns and R. Klette, “Local and global integration of discrete vector fields,” in Advances in Computer Vision, F. Solina, W. Kropatsch, R. Klette, and R. Bajcsy, eds. (Springer, 1997) pp. 149–158.

14. R. I. Mclachlan, G. R. W. Quispel, and N. Robidoux, “Geometric integration using discrete gradients,” Philos. Trans. R. Soc. Lond. A 357, 1–26 (1998).

15. S. W. Bahk, “Highly accurate wavefront reconstruction algorithms over broad spatial-frequency bandwidth,” Opt. Express 19(20), 18997–19014 (2011). [CrossRef] [PubMed]

16. X. Zhang, “An algorithm for calculating water surface elevations from surface gradient image data,” Exp. Fluids 21(1), 43–48 (1996). [CrossRef]

17. S. Ettl, J. Kaminski, M. C. Knauer, and G. Häusler, “Shape reconstruction from gradient data,” Appl. Opt. 47(12), 2091–2097 (2008). [CrossRef] [PubMed]

18. M. Grédiac, “Method for surface reconstruction from slope or curvature measurements of rectangular areas,” Appl. Opt. 36(20), 4823–4829 (1997). [CrossRef] [PubMed]

19. C. Elster, J. Gerhardt, P. Thomsenschmidt, M. Schulz, and I. Weingartner, “Reconstructing surface profiles from curvature measurements,” Optik (Stuttg.) 113(4), 154–158 (2002). [CrossRef]

20. D. L. Fried, “Least-square fitting a wave-front distortion estimate to an array of phase-difference measurements,” J. Opt. Soc. Am. 67(3), 370–375 (1977). [CrossRef]

21. R. H. Hudgin, “Wave-front reconstruction for compensated imaging,” J. Opt. Soc. Am. 67(3), 375–378 (1977). [CrossRef]

22. R. H. Hudgin, “Optimal wave-front estimation,” J. Opt. Soc. Am. 67(3), 378–382 (1977). [CrossRef]

23. B. R. Hunt, “Matrix formulation of the reconstruction of phase values from phase differences,” J. Opt. Soc. Am. 69(3), 393–399 (1979). [CrossRef]

24. J. Herrmann, “Least-squares wave front errors of minimum norm,” J. Opt. Soc. Am. 70(1), 28–35 (1980). [CrossRef]

25. W. H. Southwell, “Wave-front estimation from wave-front slope measurements,” J. Opt. Soc. Am. 70(8), 998–1006 (1980). [CrossRef]

26. M. Harker and P. O’Leary, “Least squares surface reconstruction from measured gradient fields,” in IEEE Conference on Computer Vision and Pattern Recognition (2008) pp. 1–7.

27. B. Jähne, M. Schmidt, and R. Rocholz, “Combined optical slope/height measurements of short wind waves: principle and calibration,” Meas. Sci. Technol. 16(10), 1937–1944 (2005). [CrossRef]

28. G. Strang, Introduction to Applied Mathematics, 1st ed. (Wellesley-Cambridge Press, 1986).

29. T. Davis, “Algorithm 915, SuiteSparseQR: Multifrontal multithreaded rank-revealing sparse QR factorization,” ACM T. Math. Software 38, (2011).

30. T. Davis, “SuiteSparse,” http://www.cise.ufl.edu/research/sparse/SuiteSparse/.

31. N. J. Higham, Accuracy and Stability of Numerical Algorithms (SIAM, 1996).

32. Å. Björck, Numerical Methods for Least Squares Problems (SIAM, 1996).

Method of surface topography retrieval by direct solution of sparse weighted seminormal equations

Abstract

Introduction

Mathematical formulation

Surface elevation and slopes

Arrays of slope data

Matrix-vector formulation

Rank

A priori constraints

Weighted least squares

Solution method

Synthetic data retrieval

Water topography retrieval

Conclusions

Acknowledgments

References and links:

Supplementary Material (1)

Cited By

Figures (7)

Equations (16)

Optics Express