[May 14, 2013] Water Surface Topography Retrieved from Color Images* [Journal of Atmospheric and Oceanic Technology] Tweet (Journal of Atmospheric and Oceanic Technology Via Acquire Media NewsEdge) ABSTRACT Submerged objects viewed through wavy water surfaces appear distorted by refraction. An imaging system exploiting this effect is implemented using a submerged planar light source designed so that color images reveal features of small-amplitude waves in a wind-wave tank. The system is described by a nonlinear model of image formation based on the geometry of refraction, spectral emission from the light source, radiative transfer through the water and surface, and camera spectral response. Surface normal vector components are retrieved from the color image data using an iterative solution to the nonlinear model. The surface topography is then retrieved using a linear model that combines surface normal data with a priori constraints on elevation and curvature. The high-resolution topographic data reveal small-amplitude waves spanning wavelength scales from capillary through short gravity wave regimes. The system capabilities are demonstrated in the retrieval of test surfaces, and of a case of wind-driven waves, using data collected at high spatial and temporal resolution in a wave tank. The approach of using a physical model of image formation with inverse solution methods provides an example of how surface topography can be retrieved and may be applicable to data from other similar instruments. 1. Introduction Roughness is a familiar, ever-changing feature of wind-blown water surfaces. Remote sensing techniques employed over oceans, for example, radar (Freedman et al. 2006) and lidar backscatter (Hu et al. 2008), and coupled ocean-atmosphere radiative transfer models (Chowdhary et al. 2012; Zhai et al. 2010; Ottaviani et al. 2008b; Chowdhary et al. 2006; Jin et al. 2006; Gjerstad et al. 2003) rely on models in which surface roughness is described by wind-parameterized slope distributions (Cox and Munk 1956). Detailed calculations show that for some angles of incidence and exitance, radiance reflected from rough water surfaces can have topographic dependence on features that are not captured in slope distributions (Kay et al. 2011). Of particular interest are the wind-driven waves with wavelength scales of 1 mm to 1 m, which give the most significant contribution to surface roughness (Munk 2009). Despite such interest, high-resolution topographic data are scarce, and thus there exists motivation to improve upon the methods by which topography can be reliably measured. Many instruments that have been designed to study rough water surfaces were aimed at making localized surface slope measurements. Slope measurements in a two-dimensional array of points can be retrieved from image data, for example, via polarimetric images of reflected light (Zappa et al. 2008). The use of refracted light from a submerged source to measure slope first proposed by Cox (1958) has evolved into imaging systems using measurements of light intensity, for example, monochrome images that capture slope in one dimension (Keller and Gotwols 1982), and color images for both upwind and crosswind slope vector components (Zhang and Cox 1994; Jäahne and Schultz 1992). One challenge to interpreting such images is that, while pixel response is largely sensitive to slope, it is generally a nonlinear function of both slope and elevation. Thus, interpretation of color images for slope measurement requires an estimate of the surface topography. Several techniques exist to estimate the topography from fields of slope measurements. Estimates by path integration of noisy slope data can be prone to the global propagation of error (Wei and Klette 2007). Estimates of elevation by Fourier methods can be subject to spurious high-frequency oscillations at the edges of the image (Rocholz 2008; Zhang 1996). Elevation can also be measured by including light at wavelengths that are strongly absorbed in water (such as near-infrared wavelength bands), which provide sensitivity to pathlength in addition to slope components (Jäahne et al. 2005). A pathlength and slope measurement at a point can be directly linked to the elevation, although when adjacent points are treated independently, the resulting topography can appear unnaturally rough. Presented in this paper is a system for collecting color images of a wavy surface, and a method to retrieve the surface topography from color image data. The system improves upon a setup exploited previously for the characterization of specular reflection from waves (Ottaviani et al. 2008a,c). First, a forward model of color image formation (including nonlinear effects) is used to construct synthetic images for a trial topographic solution. Pixels in the synthetic image are compared directly with pixels in a data image to provide a scalar-valued residual color error that quantifies how well each trial slope solution explains the observations.Next, the slope estimates are combined with a priori estimates of elevation and curvature to form a sparse overdetermined system of algebraic equations expressed as a linear function of topography. Least squares solutions that minimize the weighted slope, elevation, and curvature estimates are computed directly to produce an improved topographic solution (Koskulics et al. 2012). Since images have nonlinear dependence on both topography and slope, alternating slope and topography estimates are computed in an iterative process that repeats until changes in the topography become insignificant. The method is validated by measuring cases in which the topography is known: the flat surface of calm water, and a ''barge'' surface that is constructed of a thin transparent plastic sheet that is curved to draw the water surface into a sinusoidal shape. The method is applied to a case of wind-roughened surface waves in a sequence of images collected at high speed in a wind-wave tank. The topographic data are produced on a grid with 0.37-mm spacing between rows and columns extending over a 0.47 3 0.38 m area, providing the positions of over one million vertices on the surface. The resolution and spatial extent are sufficient to resolve simultaneously smallamplitude capillary and short wavelength gravity waves. The apparatus has particular limitations that are likely to restrict its usage to measurements of small amplitude, short wavelength waves in a laboratory setting. However, the analysis method provides an example of an approach that may be ideal for interpreting data from slope imaging instruments, some of which may be capable of collecting data in natural settings. 2. Experimental setup The experiments were conducted at the National Aeronautics and SpaceAdministration (NASA)Air-Sea Interaction Research Facility (NASIRF) on Wallops Island, Virginia (Long 1992). The 18.3-m interaction section of the wave tank has a 1.2 3 0.9 m (HxW) cross section as shown in Fig. 1. The water filled the tank to a depth of 0.8 m, leaving a space of air 0.4 m above the mean water level. Capabilities include air velocity adjustable from 0 to 18 m s21, and current in both directions adjustable from 0 to 5 cm s21. Air enters the interaction section from the right after passing through a diffuser to produce laminar flow at the air inlet. Beneath the air inlet is a hydraulically driven paddle whose displacement can be controlled by an arbitrary waveform generator to generate gravity waves. Located near the air outlet on the leftis a ''beach'' section to dampen waves and attenuate reflection from the tank end. The topographic imaging system, placed at the location shown in Fig. 1 approximately 8 m from the air inlet, consists of a camera mounted to a beam above the tank centerline, and a backlight with a color filter placed on the bottom of the tank. a. Camera Images are recorded by a high-speed Mikrotron EoSens camera mounted on a beam running along the tank centerline 3 m above the tank bottom. The camera is equipped with a complementary metal oxide semiconductor (CMOS) sensor covered with a Bayer pattern color filter with 1024 3 1280 pixels at 10-bit depth, and fitted with a Zeiss 85-mm lens with an f/1.4 maximum numerical aperture. Images are collected with the aperture fully open and with the exposure time set to 800 ms to minimize motion blur. The focus was set to the mean water level at a distance of 2.2 m from the camera, giving a pixel spacing of 0.37 mm and a field of view of 0.47 3 0.38 m, with the longer dimension aligned with the wind. Assuming an acceptable circle of confusion with a maximum diameter equal to the CMOS pixel spacing (14 3 14 mm), the depth of field is approximately 5 cm. In this particular configuration, the limited depth of field restricted the camera to be positioned directly above the center of the backlight so that the focal plane would coincide with the mean water level. In principle, a larger depth of field would permit alternative camera positions and enable oblique angles of viewing. The standoffdistance from the camera to the water surface could be extended or shortened using lenses of different focal length while maintaining the focal ratio. Practical limits on available lens diameters restrict the maximum standoffdistance to perhaps a few tens of meters. b. Backlight A large backlight assembly, consisting of a watertight aluminum enclosure housing 24 high-intensity fluorescent tubes and electronic ballasts and a thick glass lid covering an opening with a size of 0.99 3 0.84 m, was placed on the tank bottom beneath the camera. The glass lid is covered by a color filter consisting of plastic layers printed in cyan, magenta, yellow, and black (CMYK) ink using a commercial Hewlett-Packard large-format inkjet printer. Four layers of printed plastic were laminated together to achieve appreciable color contrast. A translucent sheet of plastic is overlaid on the color filter to diffuse the image of the fluorescent tubes so that the spectral radiance emitted from the surface is approximately isotropic. The thickness of the box is 0.17 m so that when submerged, the light-emitting surface resides 0.63 m beneath the mean water level. The color filter was produced from a CMYK pattern that was designed to optimize the color contrast of this imaging system. A test pattern, shown in Fig. 2a, was used to produce a filter having ink levels spanning the printable CMYK space in three rows and columns of perpendicularly oriented, linear gradients of cyan and yellow ink levels overlaid on nine square patches having discrete magenta levels. An image of this test pattern image was collected with the filter in place on the backlight surface while submerged in water at the mean elevation, shown in Fig. 2b. The achievable red, green, and blue (RGB) response to the full range of printed CMYK values is quite limited and results in the low color contrast appearance. Pixels having a particular level of blue response are identified in Fig. 2c and are shown rearranged on a linear red-green response space in Fig. 2d. CMYK values are modeled as a function of a red-green response using a second-order polynomial with coefficients determined by least squares regression through the selected points. To fill in the red-green space, interpolant CMYK values are found by evaluating the polynomial with the results shown in Fig. 2e. Pixels inside the limits of valid interpolation [(0#C#255) and (0#M#255) and (0#Y #255)] appear colored with all others shown in black. The highlighted area shows the largest rectangular area that can fit within the valid interpolated points and has the desired aspect ratio to fit the backlight. The particular blue response level was selected to maximize this area. The color filter is designed so that when viewed by the camera through flat water at mean elevation, images have linear gradients of red and green response, oriented perpendicularly in the upwind and crosswind directions, respectively, and a blue response that is constant across the field. The design intent is for the red and green signals to independently indicate changes of the upwind and crosswind components of the water surface slope, respectively, and for the blue signals to indicate changes in backlight output. The choice of red and green gradients provided the highest contrast images, given the limited blue response. Alternative color schemes such as polarradial gradients or hue-saturation value color schemes may also be consistent with this method, provided that the patterns are smooth and free of sharp discontinuities. The CMYK pattern selected to produce the filter is shown in Fig. 2f. Example images of calm and wavy water surfaces are shown in Fig. 3. Details on the filter performance in the context of the image formation model are discussed in section 5. 3. Image formation model The model of image formation consists of four main parts: geometric optics, backlight spectral emission, radiative transfer through the water and surface, and camera response. The geometric optics model maps the path of rays for each pixel in an image from the camera to the surface, where it is refracted according to Snell's law to a position on the backlight that depends on the local slope. The backlight emission model assumes that emitted spectral radiance is a function of position on the backlight and is isotropic in the upward direction. The radiative transfer model describes how the radiation is attenuated in the water by absorption and scattering and transmitted through the water-air boundary. The camera spectral response is calculated as an integral of the transmitted spectrum over wavelength weighted according to the CMOS sensor spectral sensitivity. Synthetic color images are constructed by repeating these calculations for each ray, corresponding to every pixel in an image. a. Geometric optics The directional transfer of light from the backlight to the camera is modeled with a single ray for each pixel. For simplicity, the model is essentially that of an idealized ''pinhole camera'' that neglects the finite aperture and diffraction effects. In principle, a more sophisticated model could be implemented to capture features such as depth-of-focus effects. To limit computation to the rays that are intercepted by the camera, the rays are traced from the camera to the backlight, opposite of the direction of radiant energy flow, as shown in Fig. 4. A Cartesian coordinate system is defined using the planar backlight surface as the reference for zero elevation, z 5 0. The origin of horizontal coordinates, x5y50, is set on the nadir axis pointing downward from the center of the camera aperture. The camera pinhole position is therefore raperture 5[0 0 zaperture]. The image rows align with the x axis, which also corresponds to the direction parallel to the wind. The camera is focused on the calm water surface with elevation z. The camera is aimed so that the nadir point can be located at some fractional row-column position [m0 n0] within the image. This position is estimated from the image of a dot placed at the point of a plumb line hung directly beneath the center of the aperture. Each pixel is highly sensitive to radiance emerging from an area of finite extent at the focal plane whose center is termed the ''pixel focal position.'' The image of a flat water surface is assumed to be undistorted so that pixel focal positions are arranged with uniform horizontal spacing Dx and Dy in the upwind and crosswind directions, respectively, estimated from the image of a ruler placed at the focal plane. Pixel focal positions of image row m and column n are rmn 5[ xmn ymn z ]5 [(m2m0)Dx (n2n0)Dy z ]. The ''pixel rays'' amn are unit vectors that point from the center of the camera aperture to the pixel focal position, written as amn 5 (rmn 2raperture)/jrmn 2raperturej. In general, pixel rays intercept a surface at an elevation and horizontal coordinates that are displaced from the pixel focal position. The horizontal displacement Dhmn is proportional to Dhmn5Dzmnhmn/(zaperture 2z), where Dzmn is the vertical displacement from the pixel focal position and hmn is the radial distance to the nadir axis .... For observations of smallamplitude waves (Dz ' 0.001 mm) at positions near the nadir axis (hmn '0.1 m) viewed from a distant camera (zaperture 2z52.5 m), this displacement is a fraction of the pixel spacing, and the horizontal coordinates are considered to be coincident with the pixel focal position. In principle, these horizontal positions could be adjusted iteratively once the change in elevation is more precisely known. Alternatively, images collected with a telecentric lens provide an orthographic viewpoint that avoids horizontal displacement entirely (Jäahne et al. 2005). The surface normal unit vector at each point of interception is computed by the cross product of vector differences between pairs of neighboring vertices, ..., as shown in the detail of Fig. 4. The direction of the refracted ray in water ... is found from Snell's law in vector form (Southall 1913): ... (1) Alternative methods for computing ... include central differences between vector pairs, ..., or the average surface normal for facets sharing a common vertex. The refracted ray intersects the backlight at some position bmn. This position can be found by computing a parameter t that scales the negative vertical component of the transmitted unit vector equal to the elevation of the surface at the point of intersection: ..., where zmn is the point of intersection and ... is the unit vector in the vertical direction, namely, .... b. Backlight spectral emission The spectral radiance emitted from the backlight was measured with an Ocean Optics USB4000 spectrometer with radiance coupled by a bare fiber probe with 100-mmdiameter aperture and 258 acceptance angle. Spectral data were measured in 3648 spectral bands spanning wavelengths from 346 to 1042 nm. A grid of uniformly spaced spectral measurements on the backlight surface in eight rows and 10 columns was collected. At the edges of the backlight, spectral radiance in all bands dropped significantly, a result of the darker edges near the end of the fluorescent tubes. After excluding rows and columns of measurements adjacent to the edges of the backlight, the radiance in each spectral band was found to vary smoothly with position. To interpolate the spectrum at points between the remaining six rows and eight columns of grid measurements, emission is modeled as a tenth-degree multivariate polynomial, including all cross terms up to the fifth degree in x and y. Inclusion of the wavelength variable in the polynomial interpolation was rejected because of the ''spiky,'' nearly discontinuous jumps of strong spectral lines thatwould require very high-degree polynomials to capture accurately. Thus, polynomial coefficients are computed by regression independently for each spectral band. The tenth-order polynomial was selected to capture small-scale variations in output while avoiding the poor numerical conditioning of higher-order regressions. c. Radiative transfer The radiance L emitted from the backlight is attenuated by absorption and scattering in water. Neglecting multiple scattering effects, the radiance transmitted from a position on the backlight to the water-air interface Lw is exponentially attenuated in wavelength l: ... (2) Here, m5z cosu is the water pathlength with u being the angle of the ray in water with respect to nadir. Spectral extinction a(l) is modeled as the sum of absorption, using the experimental data of Pope and Fry (1997), and scattering, using the model of Morel (1974). Transmittance is shown in Fig. 5a for four selected pathlengths. At the water-air interface, radiance is partially reflected, resulting in water-leaving radiance that is attenuated by a transmission coefficient T computed by the Fresnel relations: ... (3a) ... (3b) Here, Ts and Tp are the transmission coefficients for s and p polarized light, respectively; ni and nt are the refractive indices for the incident (ni 5 1.333 for water) and transmitted (nt51.000 for air) regions, respectively; and ui is the angle of incidence between ... and .... Assuming that radiance emitted from the backlight is unpolarized, T is an equally weighted sum of Ts and Tp: ... (4) For the case of barge images with plastic separating the air and water interfaces, radiance incident at angle ui upon the water-to-plastic interface is transmitted into the plastic with the coefficient of Twater-plastic(ui) given by Eqs. (3) and (4) using indices of refraction ni 5 nwater 5 1.333 and nt 5 nplastic 5 1.575. In the case of thin plastic, the layer can be approximated as having plane-parallel geometry so that the ray refracted into the plastic is incident on the plastic-to-air interface at angle uplastic 5sin21(nwater sinui/nplastic) and is transmitted with the coefficient of Tplastic-air(uplastic) given by Eqs. (3) and (4) using indices of refraction ni 5 nplastic 5 1.575 and nt5nair51.000. Themodel considered here includes the additional effects of curvature introduced by plastic with finite thickness by tracing ray paths through the plastic layer to the points of plastic-air interception where the incident angle is found. Multiple reflections between the interfaces are neglected so that the net transmission through both interfaces is given by the product: T 5Twater-plastic 3Tplastic-air. The transmission coefficient shown in Fig. 5b is nearly uniform over a broad range of ui, with maximum values of 98% from water to air and 94% from water to plastic to air. As the angle of incidence increases, transmission decreases sharply and reaches zero at 48.68, which is the critical angle of incidence for water to air. Radiance is assumed to travel from the water surface to the camera without further extinction in air. d. Camera response The color response is computed as the product of spectral radiance reaching the camera with the wavelengthdependent spectral sensitivity R(l), G(l), and B(l), for red, green, and blue pixels, respectively, estimated from a plot of data in the camera manual. Spectral radiance reaching the camera is integrated over all wavelengths forwhich the camera has appreciable sensitivity, resulting in a three-element output vector of color response for each pixel: ... (5) To increase computational speed, Eq. (5) is evaluated to form a lookup table on a uniformly spaced grid of 100 3 100 positions in x and y covering the horizontal extent of the color filter, and on a grid of 50 pathlengths m from 0.5 to 1.5 m, spanning the range of paths that are accessible by rays in this system. The grid is then used to form a lookup table that permits rapid computation of the [R G B] model output from the [x y m ] path input. Linear interpolation between grid points provides a fast method for evaluation, though grid spacing deserves attention. Coarse grids lead to model image artifacts caused by the discontinuous interpolant gradients at the boundaries between grid points. Fine grids require excessive memory. The selected grid of 100 3 100 3 50 points maintains modest memory requirements with sufficient spatial resolution so that interpolated images appear smooth. Rays having paths that fall outside of the maximum grid extent (occurring from extreme slopes that refract rays to positions outside of the backlight surface) are flagged, and the corresponding pixel color is set to a very low value near zero. 4. Topography retrieval The forward model of image formation can be approximated as a linear function for small changes in the topography about some reference state. In principle, topography estimates could be obtained directly from color data by computing linear least squares solutions using numerical solution techniques and nonlinearity can be treated by iteration. However, the large size of the resulting algebraic system, with measurement and solution vectors each having on the order of 106 elements, is computationally demanding to solve. The practicality of approaching solutions through a direct route of color to topography is limited because of the long solution times that are necessary. Direct linear least squares solutions, computed by QR factorization, is a computationally intensive step that must be repeated over several iterations for each image. Attempts of computing solutions using indirect methods, such as conjugate gradients, failed to achieve convergence despite the use of a variety of preconditioning techniques and several hours of iteration. These practical limits are overcome by dividing the problem into two simpler steps that can be solved more efficiently. First, local slopes are estimated using the color of individual pixels while holding the elevation constant using a linearized model of color as a function of local slope. This approach has the advantage of exploiting the high sensitivity of color to local slope, and it follows in the tradition of similar color imaging methods (Keller and Gotwols 1982; Jäahne and Schultz 1992; Jäahne et al. 2005; Zhang and Cox 1994; Rocholz 2008). Next, the slope estimates are combined to produce a topographic estimate. Since the coupling of slopes and topography is linear, the system remains unchanged for different slope data. The system can be decomposed into triangular matrix factors, which is an intensive computation that must be performed only once. Topographic solution vectors can be computed efficiently by forward and backward substitution. Topographic retrieval proceeds by alternating between slope and topography estimation steps, as illustrated in Fig. 6. First, slopes are estimated while holding the elevation constant. Next, the topography is estimated based on the slope estimates of the previous step. The nonlinear dependence of color on elevation is treated by subsequent slope and topography estimate iterations with the process repeating until changes in topography become insignificant. Topographic retrievals proceed by combining the slope estimates with a priori elevation and curvature estimates in a large, sparse linear system of algebraic equations. The topographic solution is computed directly by QR factorization. The procedure of slope estimation can be repeated using the improved topographic estimate to reduce forward model errors due to nonlinear elevation dependence. The procedure chosen for the cases analyzed is composed of two sequences of rotations about the x and y axes, a topographic estimate, a single sequence of x- and y-axis rotations, and a final topographic estimate. a. Surface normal vector retrieval from color Slope retrievals begin with an estimate of topography from which initial surface normal unit vectors are computed by cross product. The initial estimate can be arbitrarily chosen, such as a flat surface at mean elevation or the retrieved topography from the previous frame in a high-speed sequence, though the choice should be close to the final surface topography for rapid convergence. Slope components are retrieved by linearizing the image formation model for small rotation angles about the horizontal coordinate axes. The dotted-line box in the lower portion of Fig. 6 highlights the computational flow of slope solution for the x axis, a procedure which is then followed for the y axis. The slope estimation procedure is repeated by solving for rotations about alternating axes to treat nonlinearity. Generally, no more than two iterations are required for convergence. The small number of iterations indicates that the nonlinear model of color image formation is well approximated by linear functions. The color response for each pixel C5[R G B]T (where the superscript T indicates vector transpose) is determined by the ray path, which is a function of the pixel ray a, the surface normal vector n, and the position r of the intercepted surface, through the nonlinear model C5f (a, n, r). Vectors a and r are held constant, and the color response is treated as being approximately linear for small rotations of the surface normal unit vector from an initial reference state n0 to a new state nu 5R(du)n0, where R(du) is a rotation matrix and du is a small angle of rotation about a horizontal axis. The corresponding change dC that results from small rotations through angle du can be approximated by the linear equation dC'Jdu, where J5[ f (nu)2f (n0)]/du is a single-column Jacobian matrix. Generally, the model output of the initial reference state C0 5f (a, n0, r) will not be equal to themeasured value of color CM, resulting in a residual error vector DCM0 5CM 2C0. The linear model is inverted to determine the least squares solution of a rotation angle that minimizes the residual error: ... (6) Rotating the surface normal unit vector by this angle leads to an improved reference state, n1 5R(Du)n0, whose model color response CM1 has a smaller residual value: kDCM1k,kDCM0k. Further improvements to the slope estimate are found by repeating the rotation procedure about a different axis. For convenience, the axes of rotation are chosen as the upwind and crosswind axes. The slope retrieval procedure is performed on each pixel in an image and results in a new array of surface normal unit vectors. b. Surface topography retrieval from surface normal vectors Topography retrieval from slope estimates is accomplished by following the procedure described in detail by Koskulics et al. (2012). In this approach, the topography and slopes are related by linear finite difference equations. Topography and slopes are organized into state andmeasurement vectors, respectively, and the relationship is written as a matrix-vector product. Additional linear equations constraining elevation and curvature ... for the downwind and crosswind directions, respectively) augment this matrix. These additional equations are set equal to corresponding measurement vector elements of mean elevation and zero curvature and can be thought of as being ''virtual'' measurements. The elevation constraints ensure full system rank by eliminating the ambiguity of absolute elevation, and the curvature constraints reduce roughness caused by slope measurement error. Measurement vector elements are weighted by estimates of variance. Direct solution of the weighted seminormal equations provides an efficient means of determining the topography consistent with all available information. The topographic solution vector is computed directly. The matrix is decomposed into a product of orthogonal and upper-triangular matrices (QR factorization), with the orthogonal matrix discarded. Factoring is computationally intensive, requiring over 8 GB of random access memory and 180 s of computation time on an Intel Core i7 processor. Since the matrix remains constant, this computation is performed only once and the resulting R matrix is maintained in memory. Topographic solution vectors are computed by forward-backward substitution. This step is computationally efficient, requiring less than 5 s for each measurement vector. 5. Implementation and performance Retrieval performance is highly dependent upon the particular setup used to generate the position-dependent spectral radiance and the properties of the camera. In the setup described here, the limited color space achievable by CMYK inkjet printing processes leads to low-contrast color gradients when imaged by a camera with poor color selectivity. The dynamic range of response is further limited by the linear sensitivity range of the CMOS image sensor. Depth of field and shutter speed are other sources of image distortion. Fluctuations in backlight intensity contribute a source of systematic errors in slope retrieval. Despite such limitations, the images contain information that is sufficient to reconstruct the surface topography, as evidenced by several validation cases. a. Backlight performance Color response to radiance emerging from the backlight is a function of the spectral sensitivity of the RGB channels and radiant spectrum per Eq. (5), and can be estimated by examining sensitivity and spectral data as shown in Fig. 7. Considerable overlap in sensitivity of the three color channels results in limited spectral selectivity, which is a common feature of cameras that employ absorptive color filter elements. The red channel exhibits the highest selectivity due to the sharp drop in sensitivity between 600 and 570 nm. Thus, the red channel response is largely driven by the strong spectral emission band between 610 and 630 nm. This band is stronger on the right side of the backlight, clearly indicating an increase of red response from leftto right and resulting in the gradient of red in the direction aligned with the wind. The green channel response has an increase from the bottom to top, though the considerable overlap of green channel sensitivity with the red portion of the spectrum makes this result somewhat more difficult to discern. Cameras that use three charge-coupled device (CCD) sensors and trichroic beamsplitters provide higher selectivity. The backlight output radiance was found to change over time, though the relative spectrum remains constant. Changes in backlight output could be compensated by dividing the red and green responses by the blue response, which is largely independent of pathlength and backlight position. The two resulting quotients form a measurement vector that can convey information for both orthogonal slope components. However, the measured blue values are a factor of 2-3 smaller than red and green for two primary reasons: low sensitivity of the blue channel relative to the green and red channels, and weak short wavelength output from the backlight, both of which are evident in Fig. 7. The resulting quotients have significantly increased variance that leads to increased slope errors. The blue channel compensation was therefore rejected in this analysis, but it could be employed in improved setups. Sensitivity of the retrieved rotation angles to changes in color measurements can be estimated by examining the Jacobian matrix elements. For a pixel near the middle of the image, the Jacobian is DCM0/Dux 5[20:40 22:10 20:17 ] and DCM0/Duy 5 [22:65 20:22 20:06 ], in units of digital numbers per degree [DN (8)21], which indicates that the green and red colors channels have the highest sensitivity to rotations of the surface normal vector about the x and y axes, respectively. Inserting this Jacobian into Eq. (6) gives retrieved rotation angles that are proportional to Dux/DCM0 5[20:087 20:456 20:037 ] and Duy/DCM0 5[20:375 20:031 20:009 ](8) DN21, indicating that the green and red channels are the major determinates of rotation angles about the x and y axes, respectively. Random errors in slope are estimated from flat surface images using two iterations of slope retrievals, where the standard deviation of slopes was found to be 2.18 and 2.48 for the downwind and crosswind directions, respectively. Typical measured values of color for pixels in the center of the image are approximately 600 DN. At this level, an increase of 6 DN (a 1% change in radiance) leads to 23.58 and 22.08 changes in slope for the x and y directions, respectively. Changes on this order of magnitude, which are observed in backlight intensity fluctuations over periods of time from minutes to hours, are a source of systematic slope retrieval errors. Corrections for systematic slope errors are performed by subtracting the mean slope from each measurement prior to topographic retrieval steps. Alternative backlight intensity fluctuation correction schemes may make use of the blue channel, though this approach was not explored in great detail. b. Sources of image distortion Proper exposure settings involve a trade-offbetween exposure time and aperture f number. Long exposure time leads to motion blur caused by the time rate of change of slope, which is highest for waves with large wavenumbers and high angular velocity (small capillary waves) and where slope is discontinuous, such as the sharp crests of gravity waves. The finite exposure time made necessary by limited lighting intensity may obscure very fine features. For minimal motion blur, stroboscopic illumination is suggested.With only limited-intensity light sources of fluorescent lights, low f-number aperture settings are required. Low f-number apertures produce images subject to vignetting (reduced sensitivity) for pixels at the edge of the field of view. Low f-number settings also limit the depth of field to elevations within a few centimeters of the focal plane. Thus, sharp features on the crests and troughs of very high-amplitude waves (approximately 5 cm using this apparatus) may be blurred. The trade-offchosen here is one to minimize motion blur with the aperture fully open, and only a case of limited amplitude waves is discussed in this work. The camera attains color sensitivity with a Bayer pattern color filter mosaic overlaid on its CMOS sensor. The pattern is a rectangular array of 2 3 2 pixel unit squares with green and blue filters in the first row followed by red and green filters in the second row. Each pixel is sensitive to just one color, and three color data are obtained by ''demosaicing'' or interpolating the missing color information from neighboring pixels. For this analysis, a simple demosaicing algorithm is chosen in which the missing color data are averaged from a 3 3 3 pixel tile centered on the pixel of interest. For imageswith the pixel separation of 0.37 mm considered here, colors are sampled on a scale 3 times this length (1.1 mm). The spatially sampled image data are subject to the Nyquist criterion, which sets a minimum detectable wavelength scale of 2.2 mm. Image data may be vulnerable to distortion caused by spatial aliasing if waves shorter than the Nyquist wavelength are present with significant amplitude and spatial extent, and can be resolved by the optical system. Correctible sources of distortion include ''dark noise'' in the form of a fixed pattern, added to the radiant signal in every image. This pattern was corrected by subtracting an average of 1024 images collected with the lens cap on from the data images preceding the analysis procedure. Differences in absolute sensitivity (because of factors such as aperture f number, focal length, vignetting, exposure time, and pixel gain variance) are corrected with a multiplicative gain factor that is determined by dividing an offset-corrected flat surface data image by the model image of the surface at mean elevation. Finally, the CMOS sensor response is not linear over the full range of sensitivity, though it is approximately linear in the lower half of the range. For this reason, the exposure time was set so that the maximum response would not frequently reach beyond this threshold. However, this arrangement results in limited dynamic range and further reduces sensitivity to changes in color. 6. Validation The performance of retrieval systems can be validated by comparing retrieval results of a known standard with measurements that are made independently. Because water surfaces have dynamic and fluid properties, cases of independently measurable surfaces are considerably restricted. Examined in the following section are two cases: a calm, flat water surface and a surface created by a curved transparent plastic sheet that draws the water into a nearly sinusoidal shape (albeit with the changes due to the different refractive index of the plastic interface). As an example of real waves, a case of waves generated by a 4 m s21 wind is examined. In each of the three cases, virtual measurements are set to the mean elevation and zero curvature in both directions. Weighting is based on a 5-cm standard deviation for elevation, a very loose constraint. Experimentation with different values of curvature weighting on wavy surface retrievals showed that a standard deviation of 40 m21 provided a reasonable compromise between artificially rough retrievals and excessive smoothing that noticeably reduces small wave amplitudes. The weightings of slope measurements are chosen as the standard deviations of slopes obtained after two iterations of surface normal rotation for a flat surface case. a. Flat surface The flat surface of calm water provides a condition of known topography whose elevation should be constant, with negligible slope and curvature. Thus, ideal histograms of retrieved elevation, slope, and curvature should appear as delta functions that are centered about the mean elevation, and zero slope and curvature. Using slope data containing errors, indications of retrieval precision and accuracy can be drawn from mean and standard deviations of the retrieved results, respectively. Strong justification exists for setting the a priori variance of elevation and curvature for the flat surface to zero. However, zero variance values would lead to infinite weighting, and exceedingly small values would cause the slope measurements to have insignificant relative weighting, resulting in flat, unrevealing retrievals. To estimate the accuracy of retrievals for cases of arbitrary surfaces, for which a priori elevation and curvature variance may not be known, the weightings of all three cases are chosen to be identical. Thus, the flat surface retrieval errors can provide a rough expectation of errors for the other cases. The retrieval results for the flat surface are available in the supplemental materials (http://dx.doi.org/10.1175/ JTECH-D-12-00047.s2) and are shown in Fig. 8a. In this figure, the surface is rendered for a collimated light source using MATLAB code available in the supplemental materials (S5). The retrieved surface is free of gross defects, though some roughness is apparent. The elevation profile is shown in Fig. 8b. The histogram of retrieved elevations shown in Fig. 8c shows an irregularly shaped distribution with a standard deviation of 140 mm. Histograms of slope shown in Fig. 8d have equal standard deviations of 0.398 and 0.448 for the downwind and crosswind directions, respectively. The slope standard deviations are reduced by nearly an order of magnitude from the slope retrieval step, indicating the improvements in accuracy that can be accomplished by the least squares topographic retrieval approach. Histograms of curvature, shown in Fig. 8e, have standard deviations of 8.9 m21 for both the downwind and crosswind curvatures, again considerably less than the a priori constraints. Slope and curvature histograms appear to reproduce the bell curve shape of a Gaussian distribution. b. Plastic sinusoidal surface The sinusoidal surface consists of a thin transparent sheet of plastic curved to approximate an ideal function, in this case, a one-dimensional sinusoid with 0.5-cm amplitude and 6.5-cm wavelength. The edges of the surface are housed in a watertight frame to form a barge that floats on the water surface. When air is removed from beneath the plastic sheet, the water is drawn into the sinusoidal shape. A profile photo of this barge is shown in Fig. 9a. The RGB image used for the analysis is shown in Fig. 9b. Parts of the image outside of the sinusoidal surface, including the walls of the barge, contain pixels that are considered invalid and from which no usable information can be obtained. The area of valid data area, shown outlined in the image, is limited in length to five sinusoidal cycles that begin and end at nearly identical phase values close to the zero-elevation level. Slopes of invalid pixels are set to zero. Systematic errors in slope are reduced by subtracting the mean of the x- and y-slope data in the uncorrupted area from the uncorrupted data. The retrieval results are available for analysis in the supplemental materials (http://dx.doi.org/ 10.1175/JTECH-D-12-00047.s3) and can be rendered using a collimated light source with the MATLAB code available in the supplemental materials (http://dx.doi. org/10.1175/JTECH-D-12-00047.s5). Retrieval errors are shown in Fig. 9c. The black line indicates the location of a profile shown in Fig. 9d. The profile resembles the ideal shape (dotted line), with a close match in wavelength and phase. Elevation errors of several millimeters are apparent and are a significant fraction of the 5-mm ideal sinusoidal amplitude. Errors appearing outside the boundaries of the valid data are caused by the least squares fitting procedure, which minimizes the curvature at the boundary of valid and zero-slope (invalid) data. Retrieval errors may be partly attributed to effects that are not included in the model, such as imperfections on the plastic surface, for example, fingerprints or scratches. For this reason, it is not clear to what extent similar errors impact ordinary water surface retrievals. In principle, higher precision retrievals of arbitrarily shaped plastic surfaces may be possible with an improved imaging system. If retrieval results with higher precision could be demonstrated, despite the complexity introduced by the plastic surface, then stronger claims of validation for simpler plain water surfaces could be supported. c. Wavy surface Elevation was retrieved from a sequence of highspeed video images collected at a rate of 200 frames s21 (a video of the full sequence of renderings is included in http://dx.doi.org/10.1175/JTECH-D-12-00047.s1). The retrieval results, shown rendered in Fig. 10a using the same model as for Fig. 8a, appear similar to natural wind-driven waves. This figure is the first frame of a video sequence, available in the supplemental materials (http://dx.doi.org/10.1175/JTECH-D-12-00047.s4) and can be rendered with a collimated light source using the MATLAB code available in the supplemental materials (http://dx.doi.org/10.1175/JTECH-D-12-00047.s5). Waves of short length have significant slopes that create the contrast evident in the rendered images, but are not evident in the elevation profile shown in Fig. 10b because of their small relative amplitude. A small, darkened dimple can be seen in the bottom-right corner of the image, moving from leftto right in the video sequence. Examination of the raw data images shows a corresponding dark spot, most likely a small particle suspended in the water and drifting with the subsurface current, which was set to flow in the upwind direction. The elevation histogram in Fig. 10c has a standard deviation of 3.6 mm, indicating small amplitude waves. A pattern of positive skewness is evident, which is consistent with observations of the particular shape of water waves under the action of wind that exhibit sharp crests covering a proportionally smaller area than the broader, flatter troughs (Huang et al. 1990). The slope histograms of Fig. 10d have standard deviations of 7.58 (upwind) and 4.68 (crosswind). The smaller crosswind variance is consistent with wave crests that tend to be oriented parallel to the wind direction that are observed in Fig. 10a. The upwind slope distribution is skewed negative (surface normal vectors tilted toward the downwind direction) and is consistent with observations that these shortest waves are usually found on the downwind face of the larger waves. Crosswind slopes approximate a Gaussian distribution, while upwind slopes are highly irregular and asymmetric. Such asymmetry has long been observed (Cox and Munk 1956; Long et al. 1994; Munk 2009). Histograms of curvature in Fig. 10e also show crosswind/upwind asymmetry with standard deviations of 35 m21 (upwind) and 19 m21 (crosswind). Crosswind curvature appears symmetric about zero, though upwind curvature appears skewed negative (concave). These data, and other cases, will be examined in greater detail in future studies. 7. Summary and conclusions Color images of a submerged color planar light source refracted through a wavy water surface provide measurable signals that carry detailed information about the water surface topography. A nonlinear model of color image formation enables comparison of data images with synthetic images of trial surface topographies. A two-step iterative approach of retrieving slopes from color, followed by retrieval of elevation from slopes, enables estimation of topography from the color image data. Particular emphasis is given to the description of the color imaging system and an evaluation of its performance in the context of the image formation model. In this construction, the limited color contrast achievable by the CMYK printing process used to produce the filter and the poor color selectivity of the camera interact to limit the sensitivity to slopes. Instability in the output of the fluorescent tubes used in the backlight provides a source of systematic error in slope retrievals that is potentially correctible using normalization channels. Suggested performance improvements include alternative illumination sources, such as arrays of light-emitting diodes, pulsed light output with high-intensity stroboscopic illumination, and cameras that have better color selectivity, such as cameras with three independent CCD sensors and trichroic beamsplitters. Despite the imperfections, highly detailed topographic retrievals from color image data are demonstrated. Retrievals of a flat water surface and a sinusoidal-shaped plastic surface provide estimates of retrieval accuracy and precision. A case of topographic retrievals of winddriven waves demonstrates the efficacy of this approach in a sequence of high-speed images (see also the supplementalmaterials). Each topographic retrieval result provides a representation of the water surface in fine detail. Finally, the forward model and inverse method outlined in this paper suggest a general approach for measuring water surface topography using slope-sensitive image data based on similar or alternative techniques (e.g., slope-sensitive polarimetric or reflective color images). Such an approach has the advantage of making efficient use of all available slope and a priori data, and can account for the complexities that may be introduced by nonlinearities in the physics of image formation. While the apparatus described herein suffers from particular shortcomings, it is conceivable that an improved measurement technique could be developed and deployed at sea. The analysis methods developed here represent a step toward improved experimental methods aimed at improving our understanding of small-scale windgenerated waves on the ocean surface. Acknowledgments. The authors acknowledge support from the NASA Radiation Sciences Program and wish to thank the program manager, Dr. Hal Maring. SRL would like to express his thanks to Dr. Eric Lindstrom of the NASA Science Mission Directorate for supporting the research efforts of NASIRF during its final years of operation. * Supplemental information related to this paper is available at the Journals Online website. Supplemental information is at http://dx.doi.org/10.1175/JTECH-D-12-00047.s1 REFERENCES Chowdhary, J., B. Cairns, and L. D. Travis, 2006: Contribution of water-leaving radiances to multiangle, multispectral polarimetric observations over the open ocean: Bio-optical model results for case 1 waters. Appl. Opt., 45, 5542-5567. _____,_____, F.Waquet, K. Knobelspiesse, M. Ottaviani, J. Redemann, L. Travis, and M. Mishchenko, 2012: Sensitivity of multiangle, multispectral polarimetric remote sensing over open oceans to water-leaving radiance: Analyses of RSP data acquired during the MILAGRO campaign. Remote Sens. Environ., 118, 284-308. Cox, C., 1958: Measurement of slopes of high-frequency wind waves. J. Mar. Res., 16, 199-225. _____, and W. Munk, 1956: Slopes of the sea surface deduced from photographs of sun glitter. Bull. Scripps Inst. Oceanogr., 6, 401-488. Freedman, A., D. McWatters, and M. Spencer, 2006: The Aquarius scatterometer: An active system for measuring surface roughness for sea-surface brightness temperature correction. 2006 IEEE International Symposium on Geoscience and Remote Sensing, IEEE, 1685-1688. Gjerstad, K. I., J. J. Stamnes, B. Hamre, J. K. Lotsberg, B. Yan, and K. Stamnes, 2003: Monte Carlo and discrete-ordinate simulations of irradiances in the coupled atmosphere-ocean system. Appl. Opt., 42, 2609-2622. Hu, Y., and Coauthors, 2008: Sea surface wind speed estimation from space-based lidar measurements. Atmos. Chem. Phys., 8, 3593-3601. Huang, N. E., C. C. Tung, and S. R. Long, 1990: The probability structure of the ocean surface. The Sea, B. LeMehaute and D. M. Hanes, Eds., Ocean Engineering Science, Vol. 9, John Wiley and Sons, 335-366. Jäahne, B., and H. J. Schultz, 1992: Calibration and accuracy of optical slope measurements for short wind waves. Optics of the Air-Sea Interface: Theory and Measurement, L. Estep, Ed., International Society for Optical Engineering (SPIE Proceedings, Vol. 1749), 222-233. _____, M. Schmidt, and R. Rocholz, 2005: Combined optical slope/height measurements of short wind waves: Principle and calibration. Meas. Sci. Technol., 16, 1937, doi:10.1088/ 0957-0233/16/10/008. Jin, Z., T. P. Charlock, K. Rutledge, K. Stamnes, and Y. Wang, 2006: Analytical solution of radiative transfer in the coupled atmosphere-ocean system with a rough surface. Appl. Opt., 45, 7443-7455. Kay, S., J. Hedley, S. Lavender, A. Nimmo-Smith, 2011: Light transfer at the ocean surface modeled using high resolution sea surface realizations, Opt. Express, 19, 6493-6504. Keller, W. C., and B. L. Gotwols, 1982: Two-dimensional optical measurement of wave slope. Appl. Opt., 22, 3476-3478. Koskulics, J., S. Englehardt, S. Long, Y. Hu, and K. Stamnes, 2012: Method of surface topography retrieval by direct solution of sparse weighted seminormal equations. Opt. Express, 20, 1714- 1726. Long, S. R., 1992: NASA Wallops Flight Facility air-sea interaction research facility. NASA Reference Publ. 1277, 32 pp. _____, N. E. Huang, E. Mollo-Christensen, F. C. Jackson, and G. L. Geernaert, 1994: Directional wind wave development. Geophys. Res. Lett., 21, 2503-2506. Morel, A., 1974: Optical properties of pure water and pure sea water. Optical Aspects of Oceanography, N. G. Jerlov and E. S. Nielson, Eds., Symposium on Optical Aspects of Oceanography, Vol. 1, Academic Press, 1-24. Munk, W., 2009: An inconvenient sea truth: Spread, steepness, and skewness of surface slopes. Annu. Rev. Mar. Sci., 1, 377-415. Ottaviani, M., C. Merck, S. Long, J. Koskulics, K. Stamnes, W. Su, and W. Wiscombe, 2008a: Time-resolved polarimetry over water waves: Relating glints and surface statistics. Appl. Opt., 47, 1638-1648. _____, R. Spurr, K. Stamnes, and W. Li, 2008b: Improving the description of sunglint for accurate prediction of remotely sensed radiances. J. Quant. Spectrosc. Radiat. Transfer, 109, 2364- 2375. _____, K. Stamnes, J. Koskulics, H. Eide, S. R. Long, W. Su, and W. Wiscombe, 2008c: Light reflection from water waves: Suitable setup for a polarimetric investigation under controlled laboratory conditions. J. Atmos. Oceanic Technol., 25, 715-728. Pope, R.M., and E. S. Fry, 1997: Absorption spectrum (380-700 nm) of pure water. II. Integrating cavity measurements. Appl. Opt., 36, 8710-8723. Rocholz, R., 2008: Spatio-temporal measurement of short winddriven water waves. Ph.D. dissertation, Institute of Environmental Physics (UIP), University of Heidelberg, 164 pp. [Available online at http://hci.iwr.uni-heidelberg.de/publications/ dip/2008/Rocholz_2008_Diss.pdf.] Southall, J. P. C., 1913: The Principles and Methods of Geometrical Optics, Especially as Applied to the Theory of Optical Instruments. 2nd ed. Macmillian, 626 pp. Wei, T., and R. Klette, 2007: On depth recovery from gradient vector fields. Algorithms, Architectures and Information Systems Security, B. B. Bhattacharya et al., Eds., Statistical Science and Interdisciplinary Research, Vol. 3, World Scientific, 75-96. Zappa, C. J., M. L. Banner, H. Schultz, A. Corrada-Emmanuel, L. B. Wolff, and J. Yalcin, 2008: Retrieval of short ocean wave slope using polarimetric imaging. Meas. Sci. Technol., 19, 055503, doi:10.1088/0957-0233/19/5/055503. Zhai, P.-W., Y. Hu, J. Chowdhary, C. R. Trepte, P. L. Lucker, and D. B. Josset, 2010: A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface. J. Quant. Spectrosc. Radiat. Transfer, 111, 1025-1040. Zhang, X., 1996: An algorithm for calculating water surface elevations from surface gradient image data. Exp. Fluids, 21, 43-48. _____, and C. S. Cox, 1994:Measuring the two-dimensional structure of a wavy water surface optically: A surface gradient detector. Exp. Fluids, 17, 225-237. JEFFREY KOSKULICS AND STEVEN ENGLEHARDT Stevens Institute of Technology, Hoboken, New Jersey STEVEN LONG Air-Sea Interaction Research Facility, Wallops Flight Facility, NASA GSFC, Wallops Island, Virginia YONGXIANG HU NASA Langley Research Center, Hampton, Virginia MATTEO OTTAVIANI NASA Goddard Institute for Space Studies, New York, New York KNUT STAMNES Stevens Institute of Technology, Hoboken, New Jersey (Manuscript received 26 February 2012, in final form 2 December 2012) Corresponding author address: Jeffrey Koskulics, Light and Life Laboratory, Stevens Institute of Technology, 1 Castle Point on Hudson, Hoboken, NJ 07030. E-mail: [email protected] (c) 2013 American Meteorological Society [ Back To TMCnet.com's Homepage ]