Reconstruction of Cloud Geometry from High-Resolution Multi-Angle Images

We consider the reconstruction of the interface of compact, connected"clouds"from satellite or airborne light intensity measurements. In a two dimensional setting, the cloud is modeled by an interface, locally represented as a graph, and an outgoing radiation intensity that is consistent with a diffusion model for light propagation in the cloud. Light scattering inside the cloud and the optical internal parameters of the cloud are not modeled. The main objective is to understand what can or cannot be reconstructed in such a setting from intensity measurements in a finite (on the order of 10) number of directions along the path of a satellite. Numerical simulations illustrate the theoretical predictions.


(Communicated by John Schotland)
Abstract. We consider the reconstruction of the interface of compact, connected "clouds" from satellite or airborne light intensity measurements. In a two-dimensional setting, the cloud is modeled by an interface, locally represented as a graph, and an outgoing radiation intensity that is consistent with a diffusion model for light propagation in the cloud. Light scattering inside the cloud and the internal optical parameters of the cloud are not modeled explicitly. The main objective is to understand what can or cannot be reconstructed in such a setting from intensity measurements in a finite (on the order of 10) number of directions along the path of a satellite or an aircraft. Numerical simulations illustrate the theoretical predictions. Finally, we explore a kinematic extension of the algorithm for retrieving cloud motion (wind) along with its geometry.
1. Introduction. This paper concerns the reconstruction of a cloud surface from radiance satellite measurements. The state-of-the-art for extraction of physical properties of clouds from passive remote sensing is based on a pixel-by-pixel interpretation of the measured radiance across the solar spectrum, possibly at different angles and states of polarization. Invariably, the cloud is modeled as a plane-parallel slab [23,24]. This radical assumption about cloud geometry is more-or-less justified for stratiform clouds that have considerably more horizontal extension than thickness in the vertical. These clouds are indeed important for the balance of the Earth's climate and hydrological cycle. However, that leaves out many other important types of cloud, those resulting typically from more vigorous convection (updrafts); the cumulus class of cloud types is thus not amenable to analysis by operational remote sensing techniques. For such clouds, the first order of business in remote sensing should be to reconstruct their non-trivial outer shape, hence the goal of the present demonstration. The main objective of the paper is to reconstruct such an interface as well as the light intensity escaping from it in different directions, but without modeling the internal optical properties of the cloud.
How light propagates inside the cloud is accurately modeled by a radiative transfer model [4,18,7]. The plane-parallel slab geometry is attractive largely because it leads to one-dimensional radiative transfer, a tractable problem in computational physics [4]. As proof that one can depart radically from the plan-parallel cloud geometry and still have an analytically tractable radiative transfer problem to solve, at least in the relevant diffusion approximation, Davis [6] assumed spherical clouds; he then used his solution to derive effective optical thicknesses-or rather, diametersfor selected cumuli in a field of broken clouds observed obliquely by DOE's Multispectral Thermal Imager (MTI) [26]. Reconstruction of the constitutive (optical) parameters in a full three-dimensional radiative transfer equation is a notoriously difficult and ill-posed problem [2] that has rarely been attempted for real clouds. For computationally intense efforts in that direction during the past decade or so, see [17,25,5]. Far more efficient physics-based methods that may lead to practical implementations are currently being formulated [19,20] and demonstrated [13,14].
In this paper, we do not aim to reconstruct such parameters and rather identify geometric properties of the cloud that are directly observable from satellite measurements and are independent from its optical properties. We restrict ourselves to a two-dimensional setting to simplify the modeling and the numerical simulations. We expect most results presented here to hold without major modification in a three-dimensional environment. The cloud is modeled as a compact, connected, domain with a sufficiently smooth interface. The cloud boundary is given by a curve represented by t → γ(t) ∈ R 2 . Since the reconstructions are local, the curve is also represented locally by a graph (x, h(x)). The main objective is then to estimate h(x) from available data. To first order, the radiative transfer of light in clouds is as follows. Light is emitted from the sun, enters the cloud proper, scatters (most probably, multiple times), and then leaves the cloud at the points (x, h(x)) in all possible outgoing directions θ. Importantly for the cloud shape reconstruction, as we will show, cloud-escaping radiation is not isotropic (a.k.a. Lambertian).
The measurements obtained by sensors mounted on satellites or aircraft are modeled as the light intensity at position (x, Z) for all x and a fixed Z (the sensor platform's orbital altitude above the Earth's surface). We assume that the sensors are equipped with several narrow directional filters so that light intensity u(x, Z, θ) is measured for a discrete number of directions θ ∈ {θ j } 1≤j≤J . Typically, J = 9, as has been implemented here to mimic the Multi-Angle Imaging Spectro-Radiometer (MISR) [8].
Alexandrov et al. [1] have recently developed a completely different approach to outer cloud shape reconstruction that is adapted to a non-imaging airborne instrument, such as the Research Spectro-Polarimeter (RSP) [3]. Although not an imager, RSP is however "hyperangular" (i.e., J is several 100s). Each direction is sampled from a different point along the sensors path, resulting in a very large number of rays. In essence, all the rays that hit the edge of the cloud are used. For every pair of these rays, an ellipse is defined that fits between the paired rays, and the cloud is defined as the union of all these ellipses.
Returning to MISR, it is part of the payload on NASA's Terra satellite, launched in 1999 into a morning-crossing sun-synchronous orbit at Z = 705 km. MISR's nine imaging sensors are "push-broom" cameras that use the orbital motion pitched at fixed angles ranging from about 70.1 • from zenith (at surface) in the forward direction to the same in the backward direction. Their common spatial sampling is 0.275 km, at least in the red (672 nm) channel. Thus a typical cumulus cloud, which has commensurate horizontal and vertical dimensions of a few km, will cover a few tens of pixels in the along-track and across-track directions. In the present study, we will assume many more spatial samples. This is more like what can be achieved using a digital camera located on the ground while the cloud is advected past it by a steady wind (Z is then related to cloud height), or else by using an airborne platform (Z is then the aircraft's altitude). There are space-based imagers that can achieve very high (∼10 m) resolution, but they typically have only one view angle. Only quite recently developed airborne imaging sensors have at once high spatial resolution and multi-angle capability [9]. We will address the issues raised by the limited spatial sampling in current space-based imaging sensors in future work. In the remainder, we will continue to call the platform a "satellite" since satellites will eventually become a source of abundant free data for cloud shape reconstruction.
The outline for the rest of the paper is as follows. In Section 2, we lay out some key definitions and assumptions. The graph model for the cloud interface and the associated inverse problems are presented in Section 3. The numerical algorithm used to solve the inverse problem is given in Section 4. A more general geometric setting based on a polar representation of the cloud geometry is described in Section 5, along with numerical simulations. We summarize our findings and outline future research in Section 6. Finally, in an Appendix, we show that the available measurements are not capable of uniquely reconstructing the cloud boundary when the horizontal speed of the cloud is also an unknown.

2.
Definitions and assumptions. Let ν(x) be the outward unit normal vector at (x, h(x)), which is given with proper normalization by (1 + (h (x)) 2 ) − 1 2 times the vector (−h (x), 1). Here, h (x) is the derivative of h(x). The outward radiation at a point (x, h(x)) is therefore characterized by the intensity u(x, h(x), θ) for all vectors θ such that θ · ν(x) > 0. This is a two-dimensional function of (x, θ) whereas the available measurements correspond to J one-dimensional functions of x. This results in a severely under-determined problem and assumptions on the outgoing radiation intensity are necessary. In this study, we assume that outgoing light intensity is of the form u(x, h(x), θ) = α(x)H(θ · ν(x)) where α and H are each functions of just one variable. This physics-based model for light intensity is correct when light propagation in the cloud is modeled as a diffusive process [2]. This situation holds for sufficiently opaque clouds [7], which are also the ones for which a sharp separation between the inside and outside of the cloud is the most realistic.
That said, we are fully aware that, as hypothesized by Mandelbrot [16] and proven observationally by Lovejoy [15], cloud shapes are best represented as fractals, that is, as convoluted surfaces that do not have well-defined tangent planes nor normal vectors at any scale. This is largely due to the turbulent nature of atmospheric dynamics, especially when clouds occur. Consequently, our methodology proposes to deliver a smoothed or, better yet, appropriately-averaged approximation of the cloud's actual outer three-dimensional shape. Such smoothing is necessary in the context of limited available data, and remains a vast improvement over the aforementioned operational assumption that clouds are plane-parallel optical media, regardless of image context.
The key result of this paper is an iterative reconstruction procedure for (h, α, H), when J is sufficiently large, based on a linearization of the functional mapping (h, α, H) to the available satellite data. In "favorable" situations, which depend on the state about which the linearization is performed, this linear map is invertible and provides stable reconstructions for (h, α, H). Explicit calculations allow us to display several such favorable as well as less-favorable situations. Numerical simulations confirm the theoretical predictions.
3. Graph model for cloud reconstruction. Herein, we present the geometric assumptions on the cloud in Section 3.1. The inverse problem is then linearized in Section 3.2 and an iterative algorithm is described in Section 3.3. The calculations in Section 3.2 are the main result of this paper, as they show theoretically under which conditions some or all of the parameters of interest (h, α, H) can be reconstructed.
3.1. Geometric setting. We consider a two-dimensional bounded domain, the "cloud," with boundary given by a curve represented by t → γ(t) ∈ R 2 . To simplify the presentation, we focus here on the reconstruction of the "upper" part of the cloud, which we assume is given by the graph of a function h(x); in other words, γ(x) = (x, h(x)) ∈ R 2 , with h(x) the unknown function we wish to reconstruct. Since reconstructions will be shown to be local, an arbitrary boundary may be reconstructed by using the appropriate number of local graphs representing the interface; see Fig. 1.
The available measurements are assumed to be obtained by an imaging sensor at a fixed elevation Z; we neglect the Earth's curvature, which can be accounted for in a straightforward manner (and is in the case of actual satellites). Measurements are performed along the line X → (X, Z) and for a given number of (upward-pointing) directions θ j = (cos φ j , sin φ j ) for φ j ∈ [0, π] and 1 ≤ j ≤ J. A value of practical interest [8,9] is J = 9.
Most techniques for the reconstruction of strongly scattering optical media are based on the assumption that the equations of radiative transfer are valid to model the propagation of photons inside and outside the medium. These models require to first reconstruct spatially varying functions such as the diffusion coefficient (equivalently, the scattering mean-free-path) and the absorption coefficient. Such parameter estimations are often ill-posed. When satellite measurements are available for a large number of angles, we argue that the reconstruction of such parameters can be bypassed in some configurations, at least in a first stage. We propose an inversion procedure whose main objective is a direct reconstruction of h(x). This requires some assumptions on the solar radiation escaping from the cloud that we now present in detail.
Let u(x, z, θ) be the phase-space density of photons as a function of position (x, z) and angle θ = (cos φ, sin φ). We assume that the medium between the cloud and the satellite is clear, so that photons advect freely: This implies that for all points (X, Z) − sθ outside of the cloud, we have For a point (X, Z) on the satellite's line of measurements and for θ pointing "upwards," we denote by s(X, Z, θ) the distance from the satellite to the cloud in the Figure 1. Geometry of cloud interface direction θ and by (x(X, Z, θ), h(x(X, Z, θ))) the point of intersection on the surface of the cloud. These functions are defined by these constraints: In our theoretical analyses, we assume that the functions s and x are uniquely defined and smooth and we set s = +∞ if the half-line from (X, Z) in direction θ does not intersect the cloud. For fixed θ, we thus assume that the mapping X → x(X, Z, θ) is a local diffeomorphism.
Information from the satellite measurements thus provides knowledge of The outgoing radiation at the cloud boundary is given by u(x, h(x), θ). It is not possible to reconstruct an arbitrary radiation profile u(x, h(x), θ) from knowledge of (X, Z, θ) → u(X, Z, θ) since in fact we find a solution u(x, h(x), θ) for each choice of h(x). The outgoing radiation at the cloud surface thus needs to be constrained.
Assuming that scattering (disorder) is large in the cloud (i.e., the scattering mean free path is small compared to the overall extension of the cloud), as is often the case, then photon propagation inside the cloud is well-approximated by a steadystate 2D diffusion equation. In that case, u(x, z, θ) ∼ U (x, z) + 2θ · F (x, z) with {U, F } being the solution of equations of the form ∇·F +σ a U = 0 and F = −D∇U . This approximation breaks down at the surface of the cloud, where it needs to be replaced by the solution of a Milne problem [4,2].
At a point (x, h(x)) of the surface, let be the outward unit normal vector to the surface. Then we adopt the model with U the solution of the diffusion approximation and H an appropriate Chandrasekhar function that depends on the scattering phase function inside the cloud [4,2]. We note that Our constraint on the outgoing radiation of the cloud is to be consistent with the above-mentioned diffusion approximation. In short, we assume that for two unknown functions: α(x), a measure of the angularly integrated radiance escaping the cloud at point x; and β(θ) := H • sin(θ), an angular distribution model that is assumed independent of position in the sense that it is function only of the angle of the observed radiance with the local normal. From (3), we thus observe that with x(X, Z, θ) from (2) and ν(x(X, Z, θ)) from (4) being functionals of the unknown height function h(x). Note that multiplying α by a constant and dividing β by the same constant does not change u(X, Z, θ). We thus assume that β(φ) is normalized, for instance, by assuming that it integrates to unity. Alternatively, we can assume that its value for a given argument, such as β(π/2), is prescribed.
In summary, we thus have a measurement operator The inverse problem consists of reconstructing (α, β, h) from knowledge of M(α, β, h) = u(X, Z, θ) in (7) for a fixed value of Z, for X on a line segment, and for θ ∈ {θ j } 1≤j≤J . To illustrate, Fig. 2 shows a typical cloud model used in this study along with a set of J = 5 synthetic radiance fields based on (7).

3.2.
Linearization of the inverse problem. Let us define v := (α, β, h). The above measurement operator M, a nonlinear functional of v, does not seem to have an explicit inversion formula. Following a standard methodology, we linearize the problem about a reference configuration v 0 = (α 0 , β 0 , h 0 ) and obtain an equation for the linear update δv = (δα, δβ, δh) from knowledge of 3.2.1. Linearized system of equations. We now present the main results of the paper, namely equation (11) providing the relationship between the unknown δv and the known measurement δu.
Recall (7): u(X, Z, θ) = α x(X, Z, θ) β φ − arctan h (x(X, Z, θ)) . We wish to differentiate this expression with respect to (α, β, h). While the derivatives with respect to α and β are straightforward, the differentiation with respect to h requires a few additional steps.

Let us introduce
Differentiating the above expression for u, we find to leading order that It remains to express δx in terms of δh. From (2), we find up to second-order that We thus deduce that Then, from the above expressions for δu and δx, we obtain that The above expression is the main calculation of the paper. Our objective is to reconstruct two spatially dependent functions δα(x) and δh(x) and one angularly dependent function δβ(φ) from the spatially and angularly dependent function δu(X, Z, φ). Counting dimensions shows that the problem is formally overdetermined as soon as measurements are available for {φ j } 1≤j≤J for J ≥ 3. However, rewriting the above equation as Aδv = δu, and in normal form with A t being the adjoint operator to A, the linear operator A t A needs to be invertible. This imposes some conditions on v 0 = (α 0 , β 0 , h 0 ), which are not always met in practice. Subsequent sections present situations in which we can prove that A t A is indeed invertible provided that adapted boundary conditions for h(x) are imposed (for instance, Dirichlet conditions).

3.2.2.
Analysis of the linearized system. As recalled earlier, we assume that the map X → x 0 (X, Z, θ) for each fixed θ is a smooth change of variables. We easily verify that such is the case when h 0 (x) is sufficiently small. Then the map (X, φ) → (x 0 (X, φ), φ + ψ 0 (x 0 )) is also a smooth invertible change of variables. We then identify functions u(X, φ) and u(x 0 , φ) with φ ← φ + ψ 0 (x 0 ) and use this convenient change of notation.
Let us first consider the simplest setting with β 0 (φ) = sin φ (as a simple model for limb-darkening) and h 0 (x) = 0 (flat horizontal cloud boundary). We verify that In other words, we have If we normalize β( π 2 ) = 1, then we observe that the above expression at φ = π 2 (nadir view) provides (δ ln α)(x 0 ) or, equivalently, δα(x 0 ). It is then straightforward to observe that (δ ln β)(φ) is known up to the addition of λ cot φ for λ ∈ R arbitrary, and that (ln α 0 (x 0 )) δh − (δh) is known up to the addition of the constant λ.
Consider now the problem of the reconstruction of the parameters on a domain x min < x 0 < x max . Then all parameters are uniquely determined provided that, for instance, h(x min ) and h(x max ) are known (Dirichlet boundary conditions). Indeed, in such a setting with the additional constraint −(δh) + (ln α 0 (x 0 )) δh = S with S a known source, δh is uniquely reconstructed, and then so is δβ.
The general case with h 0 = 0 proceeds as above: we need to impose one or two conditions on δh depending on whether (1, ψ 3 , ψ 4 ) are linearly independent or not. The above analysis shows that (δα, δβ, δh) can be uniquely reconstructed with minimal additional information required on δh.
Let us conclude this section by considering the case where β(φ) is known, for instance chosen as H • sin with H the Chandrasekhar function corresponding to the appropriate scattering phase function inside the cloud. Then in the setting with cot φ and (ln β) (φ) linearly dependent, we obtain a unique reconstruction of δα and (ln α 0 (x 0 )) δh − (δh) . This provides a reconstruction of δh if for instance δh(x min ) = 0. If cot φ and (ln β) (φ) are linearly independent, then (ln α 0 ) (x 0 )δh and (δh) can both be reconstructed and δh is then uniquely determined provided that (ln α 0 ) does not vanish uniformly. Thus in practice, we expect to uniquely reconstruct (δα, δh) from the available satellite measurements.
3.3. Nonlinear inverse problem. Let us assume that the nonlinear operator A t A, as defined in the above, is invertible. Then it will be invertible for nearby values of v 0 by continuity. We may then set up the following Newton inversion for the nonlinear problem. Let us assume that v (k) = (α (k) , β (k) , h (k) ) has been reconstructed and define A (k) as the linearization of M in the vicinity of v (k) . Then up to second-order we have We thus define the iterative algorithm Provided that the initial guess v 0 is sufficiently close to v and that A t A is invertible, then it is a classical result that (A (k) ) t A (k) is also invertible and v (k) converges to v as k → ∞.
As we have seen in the preceding section, the matrix A t A may not be invertible or may be ill conditioned; for instance when h 0 (x 0 ) = 0 and α 0 (x 0 ) = 0, which imposes that ψ 3 = 0 in (13). To remedy this issue, we add a small penalty term and solve instead with B a linear operator that penalizes the curvature of h, for instance, B can be a discretization of h (x); see [10] for references on regularization. We then choose λ small in order not to affect the reconstruction of h when the curvature does not vanish.
As a final remark, we note that nonlinear problems may be injective even when their linearizations are not (think of the map x → x 3 whose linearization is not invertible at 0). However, when h(x) = h 0 and α(x) = α 0 are constant, then the available measurements M are independent of h 0 so that both the linear and the nonlinear inverse problems are not injective. 4. Numerical discretization of the graph formulation.

4.1.
Description of the algorithm. The upper boundary of the cloud can indeed be represented as a graph (x, h(x)) for a large class of clouds. However, the lateral and bottom boundaries would require different parameterizations. In order to avoid technical difficulties, we have considered clouds with the following (simplified and somewhat unrealistic) structure. We assume that the bottom of the cloud is flat and given by the line segment between (x L , h B ) and (x R , h B ). The left part of the cloud is assumed vertical and given by the line segment between (x L , h B ) and (x L , h 1 ). The right part of the cloud is also assumed vertical and given by the line segment between (x R , h B ) and (x R , h N ). The upper part of the cloud is represented by (x, h(x)) for x L < x < x R and discretized as a continuous piecewise linear function with nodes (x j , h j ) for 1 ≤ j ≤ N such that The outgoing radiant flux (angularly-integrated radiance) at each of the N − 1 linear pieces of the upper surface is denoted by α j for 1 ≤ j ≤ N − 1. The fluxes on the left and right sides are denoted by α L and α R , respectively. The function β is discretized by a continuous piecewise linear function equal to β p at the angles π p P for some given value of P equal to or comparable to J. With this geometry for the cloud, it remains to construct a discretization of the forward map M. We assume measurements available for {φ j } 1≤j≤J and discretize the upper surface by X n = n∆X for a discretization step ∆X and n ∈ Z. The value u nj is then the integral of the real solution u(X, Z, φ j ) generated by the discretized cloud averaged over the spatial interval [n∆X, (n + 1)∆X]. We verify that u nj is non-zero for a finite number of values n.
The unknown parameters and (discretized) functions are therefore while the information is given by u nj . The reconstruction of x L and x R is performed as follows. We assume that φ j = π 2 for some j and look at the support of u nj . The smallest value n L such that u n L j = 0 and the largest value n R such that u n R j = 0 define x L = n L ∆X and x R = n R ∆X. The above procedure defines a continuous functional from (h B , h n , α n , α L , α R , β p ) to u nj . That functional has been linearized as explained in the derivation of (11) and the corresponding nonlinear problem inverted as explained in Section 3.3.

Blocking of light.
For geometries with rapidly varying h(x), some photons emitted from a point x 0 in a direction θ are such that the half line 0 < s → x 0 + sθ intersects the boundary of the cloud before reaching the satellite line (X ∈ R, Z). In the numerical simulations, each half line 0 < s → x 0 + sθ j that starts from point x 0 is tested. We calculate the signed distances (d 1 , d 2 ) from every pair of consecutive points (x 1 , x 2 ) on the cloud, to the half line 0 < s → x 0 + sθ j . The intersection happens when d 1 d 2 < 0. When it intersects the cloud before reaching the satellite line, then the effect of the radiation emitted from (x 0 , θ j ) on the measurement u(X, Z) is disregarded. The linear analysis of the preceding section shows that the reconstruction of (α, β, h) is possible so long as a sufficiently large number of directions θ j reach the satellite line. This is confirmed by the numerical simulations presented below.

Examples of numerical inversions.
We now present some numerical simulations that display the performance of the algorithm. All reconstructions perform the reconstruction of α, β, and the height h. In all algorithms, β is given by the values presented in Fig. 3 and represents a sine function. The initial guess and a typical reconstruction for β are also shown on that figure.
For the rest of the section, we focus on the display of the reconstructions of α and h. In all simulations, the number of angular measurements is J = 9. The spatial discretization is given by N = 51 grid points.
In our examples, we used not only arbitrary functions for α, such as a constant or a step between two non-zero values, but also those that have more resemblance to real-world situations. Such cloud shape-dependent α functions are compiled as follows. Let −ξ be direction above the horizon (−ξ z > 0) from which the sun sheds light on the cloud. Let ν(x) be the unit normal vector on the surface of the cloud that was defined earlier, and set 0 < ρ < 1. Then, (16) α(x) = −ξ · ν(x), if −ξ · ν(x) ≥ ρ, and sun is not blocked; ρ, otherwise.
Thus {x, such that α(x) = ρ} represents the self-shaded portions of the cloud's surface, where light emerges as a diffuse field. By the same token, {x, such that α(x) ≥ ρ} represents portions of the cloud's surface directly exposed to the sun.
In the figures below, we always use the same three colors (blue, red and black) to represent the true value, the reconstructed value and the initial guess respectively. In all cases, the figures are displayed with h(x) on the left, using the same unit of length as on the horizontal axis for geometrical correctness, and α(x) on the right.
Example #1. In this case, α is chosen to be a constant, and the cloud has a moderate amount of curvature: Example #2. α is now chosen to be a step function, with the same cloud as in #1: Example #3. In this case, α is from (16) with ξ = π 2 , meaning overhead sun, and ρ = 0.2, for the same moderately curvy cloud: Example #4. Same α as in #3, but the cloud has a lot more curvature: Example #5. Same cloud as in #4, but now ξ = π 6 , meaning sun at 60 • from zenith (as indicated): In the five examples we presented above, the construction algorithm works well in situations where h(x) and α(x) are reasonably smooth, even for an α(x) that has a jump discontinuity.

5.
Reconstruction of more general geometries.

5.1.
Cloud boundaries in polar coordinates. The constructions presented so far assume that the cloud surface is represented by a graph (x, h(x)). Clouds with a closed boundary thus cannot be reconstructed by such a method unless several charts are considered, and that there are direct observations of each one. Several modifications of the theory allow us to reconstruct clouds with a closed boundary. The simplest method is arguably to assume that the boundary can be represented as a graph in polar coordinates [0, 2π) θ → r(θ) ∈ R 2 . The measurements are now performed on a circle of radius R that shares the same origin with the cloud, where R > max θ [r(θ)] is fixed and known. Again, we assume that we receive measurements for {φ j } 1≤j≤J with J ≥ 3.
For an Earth science application, one can envision here an aircraft carrying an imager that circumnavigates an opaque cloud of interest; this "cloud" could equally well be an opaque particulate plume emanating from a powerful source (e.g., ash from an erupting volcano, smoke from a large wildfire). Alternatively, one can envisage a network of ground-based cameras that observe low clouds advected by the prevailing wind. Let u(θ, r(θ), φ) be the light intensity as a function of position (θ, r(θ)) and angle φ. Assuming the same conditions as in previous section, we have where (Θ(θ, r(θ), φ), R) is the position on the circle where we receive the measurements.
Let t(θ, r(θ)) be the slope with respect to the X-axis of the tangent line at a point (θ, r(θ)) of the surface, and let (17) ν(θ, r(θ)) = 1 be the outward unit normal vector to the surface in Cartesian coordinates. Then, as before, we assume that for two more unknown functions α(θ) and β(φ). From the above, we thus observe that where we continue to assume that β(φ) is somehow normalized.
In this setting, we thus have a measurement operator The inverse problem consists of reconstructing (α, β, r) from knowledge of u(Θ, R, φ).

5.2.
Description of the algorithm. Let us represent the cloud boundary by a graph (θ, r(θ)) for θ ∈ [0, 2π), which we assume to be piecewise linear. For 1 ≤ n ≤ N + 1, we have r(θ 1 ) = r(θ N +1 ) and θ n = θ 1 + 2π n−1 N , hence ∆Θ = 2π/N . The radiance at each of the N linear pieces of the cloud is denoted by α n for 1 ≤ n ≤ N . The function β is discretized in the same way as previously.
With this geometry for the cloud, it remains to construct a discretization of the forward map K in (18). We assume measurements available for {φ j } 1≤j≤J and discretize the measurement circle of radius R by Θ n = n∆Θ for a discretization step ∆Θ and number n ∈ Z. Note that n depends on R because, the bigger the radius, the finer we have to discretize (increase n) in order to get reasonably smooth measurements. The value u nj is then the integral of the real solution u(Θ, R, φ j ) generated by the discretized cloud averaged over the arc (Θ n R, Θ n+1 R). We verify that u nj is non-zero for a non-vanishing number of values n and j.
The above procedure defines a continuous functional that maps (θ m , r m , α m , β p ) to u nj . The functional can be linearized similarly as explained in the derivation of (11) and the solution of the corresponding nonlinear inverse problem is similar to what has been explained in Section 3.3.

Numerical simulations.
We now show some numerical simulations generated from the model we just described. In our settings, we chose (0, 0) to be the origin of our cloud and of the circle of radius R where we receive the data. In our examples, we used not only arbitrary functions for α, but also the more realistic one in (16) that mimics solar illumination. In the following figures, we continue to use three colors (blue, red and black) to represent the truth, the reconstruction and the initial guess, respectively. Our initial guess for α is always a constant function α(θ) = 1. In all cases, the figures are displayed with r(θ) on the left and α(θ) on the right.
In this example, we obtain excellent reconstructions of both r(θ) and α(θ).
In this example, the reconstructions for both r(θ) and α(θ) are good except in the vicinity of the two points highlighted in green and black boxes in both figures (upper middle to the right). Note that we have mentioned earlier that the reconstruction would break down when α = 0 and h = 0, where (x, h(x)) is a Cartesian representation of the cloud surface in an appropriate system of coordinates. We could see clearly that some point between the highlighted points is precisely where α = 0 and h = 0. The reconstruction would have been worse had we not added the small penalty term described in §3.3.
In this configuration, all hypotheses necessary for a good reconstruction are met, leading to a well-reconstructed profile.
Example #4. In this case, α is chosen to be a step function.
In this last example, we also gave more curvature to the cloud (not illustrated) and again the reconstructions were good for the most part, but not in the vicinities of the points highlighted in green and black boxes, which were highlighted in both figures. We see that the reconstruction broke down for the same reason as in Example #2. 6. Summary and outlook. This paper presents reconstruction procedures for the geometry of an isolated cloud from airborne or satellite observations at J ≥ 3 viewing angles and at quite high spatial resolution (on the order of 100s of pixels across the cloud). Rather than reconstructing the optical properties of cloud, which may vary in space, we focus here on a robust reconstruction of the geometry of the cloud without modeling light scattering.
In the simplified setting of a two-dimensional cloud, we showed that the three onedimensional functions-(h, α, β), see forward model in Eqs. (6)-(8)-could all be reconstructed from multi-angle remote sensing measurements in favorable situations where the linearized map is invertible, and when J is sufficiently large. This was confirmed by several numerical simulations showing that our iterative procedure converged in many settings. However, there are situations where such a map is not invertible. Specifically, convergence problems arise when the cloud's boundary lacks curvature and/or the angular distribution of light emerging from the cloud is nearly isotropic.
More surprising, we found that the horizontal motion of the cloud, when unknown, could not always be uniquely recovered from the available measurements (see appendix). This finding immediately motivates further research since cloudbased wind speeds are routinely retrieved from multi-angle imaging sensors in space with far lower spatial resolution than assumed here (on the order of 10s of pixels across the cloud, or less) [12]. On the other hand, no attempt is made yet with these sensors to retrieve cloud shape, only their top height in some spatially-averaged sense [21].
In the present study, the synthetic measurements were produced with the same forward model as used in the reconstruction ("inverse crime"), which is acceptable in a demonstration. We plan to repeat the reconstruction using data from a more realistic forward model involving 2D (and eventually 3D) radiative transfer inside the cloud. Similarly, more realistic clouds can by used to generate the measurements. It will be particularly interesting to see how the present algorithm, which assumes a smooth cloud boundary, reconstructs a smoothed version of a fractal cloud boundary.
The main result of the paper is to show that geometric features of the cloud could be reconstructed from satellite measurements without having to estimate the radiative transfer parameters inside a cloud. In some configurations, this procedure may have too strong limitations. In such a case, a parameterization of the internal optical properties of the cloud becomes necessary. The resulting inverse problem for (h, α, β) and such new properties is the object of current research.
Although the present demonstration is in two-dimensional space, we expect many results to hold in a three-dimensional environment where we can apply the new cloud boundary reconstruction algorithm to real remote sensing data. The main new feature is that the now three-dimensional normal vector ν(x) is no longer in the plane given by the satellite trajectory and measured directions θ j , which would slightly modify the form of the linearized operator mapping δv to δu. This straightforward extension will open the way to processing real-world cloud imagery from airborne [9] or satellite [8] multi-angle sensors.
Appendix: Determination of cloud speed. Let us now assume that the cloud moves in the x direction at constant speed v c . We also assume that the satellite moves with speed v s along the X axis and that v s > |v c |. (Both velocities are assumed to be far less than the speed of light.) When v c = 0, measurements are performed at X = v s t for all values of X or, equivalently, of t. This is the setting of the preceding section. When v c = 0, then the cloud is immobile in the reference vs X. We define λ = vs−vc vs an unknown parameter when v c is unknown (we assume v s known).
The inverse problem now consists of reconstructing (α, β, h, λ) from knowledge of u. Note that λ is an additional scalar coefficient.
Then as in the preceding section, we observe that δα(x 0 ) + α 0 (x 0 )Xδλ can be reconstructed under appropriate assumptions on δh. Without prior information on δα, it is impossible to separate δα from δλ.
In the simplified case where h 0 = 0 and h 0 = 0, and considering only the part of (A.2) involving δα and δλ, we find δλ .
The functions (1, tan φ tan φ−h 0 (x0) ) are then linearly independent as soon as h 0 (x 0 ) = 0. In such a setting, δλ and δα can be separately reconstructed, as well as δh and δβ under conditions on δh similar to those obtained in the main text.
This shows that, in spite of what looks like redundant measurements (u(X, φ) known for all X and φ), the reconstruction of the additional scalar coefficient δλ is not guaranteed in general. The nonlinear problem for (α, β, h, λ) is then handled as described earlier. That all coefficients in (α, β, h, λ) cannot be reconstructed has been confirmed in numerical experiments, where we were not able to obtain any converging algorithm. The simulations (not presented here) were carried out for values of λ in the range 0.7 to 0.9.
Simultaneous reconstruction of the cloud's shape and speed requires a different modeling of the optical parameters of the cloud. This is the subject of ongoing research. To motivate this future work, we note that horizontal cloud velocity is successfully retrieved at the same time as cloud (top) height using multi-angle imagery with 275 m resolution from the previously mentioned MISR instrument on NASA's Terra satellite [11,12]. In that case, λ is very nearly unity since v s is an orbital velocity (∼7 km/s) while v c is generally only a few m/s. Apart from high-accuracy geolocation, the main requirement is to locate the same cloud or cloud feature as seen by different MISR cameras using standard feature matching methods [22], as is done routinely for cloud height estimation by stereography [21]. A minimum of three cameras is required to unravel cloud height and the along-track component of the wind.