STOCKPILE MONITORING USING LINEAR SHAPE-FROM-SHADING ON PLANETSCOPE IMAGERY

The storage and management of stockpiles of materials is a fundamental process in large scale activities such as mining, civil engineering, and in the management of waste landfill sites. Following the evolution of stocks has always been important, and advancements in remote sensing technology are not only facilitating this, but also making it possible in near real-time. Nowadays, this monitoring appears to be performed almost exclusively using UAV based techniques. This paper proposes to apply a simple Shape-from-Shading method on low resolution satellite images provided by PlanetScope to monitor the evolution of the volume of stockpiles. The proposed Shape-from-Shading formulation makes it possible to handle occluding objects in the scene. The loss in accuracy due to the low resolution of PlanetScope images is well compensated by the daily revisit frequency and by the fact that spaceborne acquisitions require no human supervision. With satellites it is also easy to follow simultaneously several stocking sites all over the world. We test our method on two coal storage sites and demonstrate that the stockpiles are well detected.


INTRODUCTION
The storage and management of stockpiles of materials is a fundamental process in large scale activities such as civil engineering, mining and extraction, and in the management of waste landfill sites. The ability to follow in real time the evolution of stocks is more and more important today. For instance, following the evolution of coal stocks over two years in European ports allowed to monitor the switch from coal to gas as countries seek to reduce their carbon footprint to combat climate change (Hodges, 2019). On a more mundane level, knowing the evolution of coal stocks on a day-to-day basis makes it possible to anticipate the evolution of prices.
All these techniques are more precise, faster, safer and cheaper than the manual surveying measurements undertaken until a few years ago. The sheer number of commercial actors proposing UAV-based solutions proves the popularity of these methods. However, they fall short when the goal is to survey a large number of stocking areas around the world with a high frequency.
Satellite images have the advantage over UAV-based imagery that they require no supervision and can have a high revisit time. They also allow surveillance in areas of difficult access. As with UAV imagery it is possible to acquire stereo pairs/triplets from high resolution satellites to recover 3D surface models (Facciolo et al., 2017, de Franchis et al., 2014. These techniques could be used to monitor the evolution of stockpiles with satisfying results. High resolution satellite stereo imagery is however expensive, which limits the applicability of this approach for continuous monitoring. New micro-satellite constellations enable unprecedented systematic monitoring applications thanks to their wide coverage and short revisit capabilities. The PlanetScope constellation (Planet Team, 2017), for instance, provides daily global coverage at a resolution of 3.5m. These types of products are not adapted to do stereo reconstruction, nevertheless the stockpiles are clearly visible in the images.
As the goal of this work is to study the evolution of stocks, micro-satellites provide a good trade-off between accuracy and revisit time. Thus, in this paper we propose to use Shape-from-Shading techniques (Horn, Brooks, 1986) on these lower resolution Planetscope images. This provides a much cheaper monitoring solution. Although the recovered shapes are less accurate than the ones acquired by stereo, the information they provide is enough to study the overall dynamics of the stocking sites (as illustrated in Fig.1). The use of mixed boundary conditions, as opposed to pure Dirichlet conditions, allows us to handle parts of the image that do not comply with the Lambertian model. We demonstrate on two sites that this technique allows us to automatically monitor stockpiles. We also include a validation of our volume computation using stereo-based methods on some of the dates.
The rest of the paper is organised as follows. In Section 2 we review the problem and the literature on shape-from-shading.
In Section 3 we present the proposed method. Experiments are shown in Section 4 and we conclude in Section 5

SHAPE FROM SHADING
The Shape-from-Shading (SFS) problem consists in recovering the three-dimensional shape of a surface from one image of that surface. One of the main interests of this problem is that the input data are minimal as a single image is used. This is an important difference with other three-dimensional reconstruction methods such as stereo.
This problem was first formulated by Horn (Horn, 1970) as finding the solution of a nonlinear first-order Partial Differential Equation (PDE) called the brightness equation Under the assumption of orthographic projection, the visible part of a scene is, up to a scale factor, a graph z = h(x, y) where (x, y) is an image point. The unit normal n can be easily expressed as The light source is assumed to be at infinity with a unit direction vector d = (α, β, γ). For a Lambertian surface of uniform albedo A, R(n(x, y)) = Ad · n(x, y).
Assuming an albedo equal to 1, (1) can be written This equation can be generalised to consider perspective projections (Lee, Kuo, 1994), a light source no longer at infinity but at the centre of projection (Prados, Faugeras, 2005) or non Lambertian scene (Bakshi, Yang, 1994). It has been proved in (Belhumeur et al., 1999) that when the lighting direction and the albedo of the surface are unknown, the same image can be obtained by a continuous family of surfaces.
A large number of numerical methods has been proposed for solving the Shape-from-Shading problem since its introduction in the 70s. Two review papers (Zhang et al., 1999, Durou et al., 2008 have extensively described these methods. They have classified them into three main classes: the methods of resolution of PDEs, the optimisation methods and the methods approximating the brightness equation. The first class contains characteristic strips expansion and approximation of viscosity solutions. The use of photogrammetric cues for 3D reconstruction was first proposed in (Van Diggelen, 1951) and was later developed in (Horn, 1975). This method has been essentially used for the theoretical study of the number of solutions of class C 2 of the eikonal equation that corresponds to a frontal light source at infinity. The notion of "viscosity solutions" is a more recent approach introduced in (Rouy, Tourin, 1992). Several algorithms have been proposed to compute these solutions, among which finite difference numerical methods (Rouy, Tourin, 1992, Lions et al., 1993 and Markov Chain approximation (Oliensis, Dupuis, 1993).
The optimisation methods are based on the variational formulations. They differ in the choice of the unknown, the regularisation term of the functional and the method of minimisation. A major reference for this class of methods is (Horn, Brooks, 1986). The method we use in our experiments belongs to this class.
Finally the third class covers the local methods and the linear methods. Local methods compute the normal at each point of the image independently of the other points. They usually assume the surface to be locally spherical (Hayakawa et al., 1994) or cylindrical (Wildey, 1986). Linear methods make a global (Kozera, Klette, 2000) or a local (Tsai, Shah, 1994) linear approximation of the reflectance function (3). A more detailed description of these different classes of methods can be found in the surveys previously mentioned.

OUR METHOD
In order to apply our method we first remove the useless images containing clouds or that are too blurry (Anger et al., 2019). Detecting clouds in images from pushbroom satellites can be achieved by exploiting the parallax effect on the different channels (Dagobert et al., 2019a). However, this method is not adapted to our input images as PlanetScope satellites do not use pushbroom sensors. But, since we are in possession of a time series of the same scene it is possible to use the algorithm proposed in (Dagobert et al., 2019b). Its approach exploits the redundancy of information acquired during revisits. On the other hand, the input images are not perfectly aligned. We therefore align them with the phase correlation method described in (Leprince et al., 2007) and studied in detail in (Rais, 2016). Once the images are aligned, we delineate the coal storage areas. This has to be done manually but only once for each site.

Continuous Shape-from-Shading Formulation
Like most shape-from-shading methods we assume a Lambertian scene, an orthographic view and a far light source. The metadata of the satellite image provide the azimuth φ and elevation θ of the sun from which the light direction d = (α, β, γ) can be computed: ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-2-2020, 2020 XXIV ISPRS Congress (2020 edition) The goal is to find a digital elevation model h(x, y) compatible with the observed data. As data, we have the orthoimage I(x, y), a domain of interest Ω and a boundary condition for h given on the boundary ∂Ω of the domain. For the sake of simplicity-and because it suits our concrete applicationwe will always assume homogeneous Dirichlet boundary conditions (h = 0 on ∂Ω); however, the techniques described below can be easily adapted to a boundary condition with arbitrary values. Furthermore, we also allow a Neumann-type boundary condition on a small part Γ Ω of the boundary, that is ∇h · n = 0, where n is the normal vector to ∂Ω.
The shape-from-shading equation for a Lambertian model with albedo A and constant ambient light B has the following form This is a nonlinear first-order equation that can be studied by the method of characteristics. In the general case it has no unique solution, or it may even fail to have a solution at all, and further conditions are required. Here we propose a particular simplification of this equation that is suitable for our application.

Linearised Shape-from-Shading
The proposed method relies on the following two assumptions: 1. The sun is low in the horizon (strictly lower than 45 • ); 2. The entire terrain is lit by the sun (no shadows).
Notice that these assumptions taken together imply that the maximal slope of the terrain is 45 • . These assumptions are of course not valid in general, but they are quite realistic for the application to heliosynchronous images of loose material stockpiles (whose angle of repose is always smaller than 45 • , attained for wet sand). These assumptions allow for a great simplification of our equation, leading to the linearised solver. The first consequence of the assumptions is that the typical "inversion ambiguity" of shape-from-shading models is not possible: you cannot invert the whole terrain profile of a solution and obtain another solution like when sun is very high. The second consequence, regarding small slopes, is that we can assume that ∇h 1 so that 1 + h 2 x + h 2 y ≈ 1, which turns the nonlinear equation into a, mostly trivial, linear equation. The proposed technique, described below for the discrete case, consists in solving this linear equation and then iteratively refine the solution by solving a sequence of analogous linear problems that approach the non-linear one.
Setting 1 + h 2 x + h 2 y = 1, the shape from shading equation (6) becomes linear. After arranging the terms we get the following problem: where α = (α, β) is the projected direction of the sun and f = ( The output of the shape-from-shading is robust to changes to the mean of f . The mean normalised image is shown in (a) and the resulting surface in (b). In (c) and (d) we modified this mean to observe that it results in a slight slope on the reconstructed surface. Although some differences can be observed along the boundaries of the bands the stockpiles are barely impacted.
The shape-from-shading equation is based on the assumption that the intensity depends only on the surface normal, illuminance and albedo, with the albedo often put to one when unknown. For our input images we can see that the intensity is also strongly dependent on the atmospheric conditions as seen in the top left columns of Figures 4 and 5.
Note that for any point on an horizontal surface we have |∇h| = 0, thus f (I, h, γ) = 0. By assuming the ground to be horizontal and the stockpiles to be more or less regular, enforcing a zero mean for f ensures the integrability condition. Figure 2 shows that the reconstruction is robust to changes of the mean. If the values corresponding to the ground in f are zero the ground appears flat in the output image. If this doesn't hold true the SFS will reconstruct the ground as a slope. This is visible only when comparing the values along the sides of the lanes of stored coal.
While we do not give any proof of the convergence or soundness of this simplified procedure for larger slopes (when the approximation ∇h < < 1 does not hold), the experiments show that it consistently produces surfaces of a reasonable shape, and quantitatively comparable to those obtained by stereo methods.

Study of the one-dimensional case
The linearised equation (7) is a first-order linear equation on a compact domain. In general, these equations never have a solu-tion unless the data satisfies a strict integrability condition. As we explain below, the integrability condition can be enforced by applying an affine contrast change to the input image. Any first-order PDE can be rewritten as a family of independent ODE along curves that cover the domain, called characteristic curves of the equation. In our case, the characteristic curves are straight line segments in the direction of the sun. Along each of these curves, the ODE to solve has the following form (we consider first the homogeneous boundary condition): where f (x) is the intensity at the position x and u(x) is the height profile to be found. In general, this equation has no solution that satisfies both boundary conditions, as can be seen by integrating on the whole domain The general case (in 1D) is thus given by an equation of the following form where R and S are constants that we can chose so that the integrability condition holds. The equation is directly integrable and the solution is Notice that u(a) = p is automatically satisfied, and the condition u(b) = q is written as q = R b a f + S(b − a), which determines the unique possible value of B: The value of R = 0 is thus arbitrary, as it amounts to fixing the (unknown) albedo of the surface. We can set it to 1 for simplicity.
To solve the linearized problem we use a variational approach with the regularisation term proposed in (Brooks, Horn, 1985). This leads to minimising an energy of the following form The Euler-Lagrange equation of this quadratic functional is The term D 2 (h)( α, α) is the second-derivative of h in the direction α, which does not "see" variations of h perpendicular to α. These variations are smoothed-out by the laplacian term, which is an isotropic regularizer. This Euler-Lagrange PDE is a linear elliptic equation that has a unique solution provided there is at least one Dirichlet datum on each connected component of ∂Ω.

The effect of occluding objects and the necessity of mixed boundary conditions
An important hypothesis of our SFS method is that the albedo is constant and that there are no shadows. This is not true in the presence of occluding objects. In our experiments these occluding objects are cranes, that tend to be much brighter than the coal stockpiles we observe. A simple thresholding of the intensity images detects them well in most cases. Once detected, their boundaries are treated differently from the manually drawn boundaries of the lanes. Assuming that the ground on which the coal is deposited is flat it is coherent to enforce Dirichlet boundary conditions. On the other hand when a crane is above a stockpile the reconstructed height should not be imposed to be zero on that part of the boundary (Dirichlet boundary conditions), but instead to have zero normal derivative (Neumann boundary conditions).
The importance of using these mixed boundary conditions is noticeable in Figure 3, particularly on the middle lane. The top row corresponds to the normalised image and the reconstruction without removing the cranes. This leads to a depression on the ground of the second lane and a peak on the bottom lane. Imposing boundary Dirichlet conditions everywhere (bottom left of the figure), the former aberrations disappear but the left stockpile appears carved in the middle. This comes from enforcing a zero height near the middle of the pile. Adding Neumann conditions on the occluded part (bottom right) leave the reconstruction more freedom along the crane boundaries and the results are much more realistic (bottom right).

Discrete formulation
To ease the discretisation of equation (7), and especially the Neumann boundary condition, we use the algebraic graph formalism for finite difference schemes (Lézoray, Grady, 2012). Briefly, the whole image domain is modelled as a grid graph of size W × H. This graph has n = W × H vertices and m = 2W H − W − H edges. Notice that m ≈ 2n. The region of interest Ω is thus a subset of the vertex set, identified by a binary vector 1Ω ∈ R n , and the boundary condition is specified by the knowledge of h ∈ R n outside of Ω. Scalar fields correspond to real-valued functions defined over the vertices (i.e. vectors of R n ), and vector fields to functions defined on the edges (i.e. vectors of R m ). We define the following three matrices: The identity matrix I ∈ Mn,n(R), the mask operator of the region of interest M = diag(1Ω) ∈ Mn,n(R), and the signed incidence matrix of the graph, ∇ = Mm,n(R). The matrix ∇ is the operator for computing first-order derivatives by forward differences.
Notice that if u ∈ R n represents a scalar field, then its discrete gradient is ∇u ∈ R m , which is a vector field according to the above convention. From these operators we define the divergence operator −∇ ∈ Mn,m(R) and the Laplacian ∆ = −∇ ∇ ∈ Mn,n(R). We also define a centering operator C = 1 2 |∇| ∈ Mn,m. It can be checked readily that the pointwise scalar product of two discrete vector fields F, G ∈ R m corresponds naturally to F · G = C (F G), where is the component-wise product (Hadamard product) on R m . This allows to compute the directional derivatives that appear in the continuous equation. Notice that the Hadamard product can be written as a matrix-vector product F G = diag(F )G. To enforce Neumann boundary conditions, we remove the lines of the matrix ∇ that corresponds to the edges that we want to remove.
This formalism allows a direct translation of the linear equations and variational problems defined above in a continuous setting. For example, the linear shape-from-shading equation f = α · ∇u is discretised as f = C diag(Sα)∇h, where Sα ∈ R m is a (constant) vector field in the direction α. This is a square linear system with matrix Nα := C diag(Sα)∇ ∈ Mn,n(R). The matrix Nα computes the discrete directional derivative in the direction of α. The Dirichlet boundary conditions g are linear constraints that can be easily added to the problem: or (M A+I −M )h = M f +(I −M )g. This is a linear problem on the unknown h. As happens in the continuous setting, this problem has no solution (its matrix is singular) unless the data term happens to satisfy a further integrability condition. Worse than that, even when there is a solution, the computations along each characteristic line are performed independently, resulting in a structured noise on the solution. To solve this problem we prefer to use the variational formulation, that corresponds to the associated normal equations of the linear problem, and where it is easier to add a regularisation term.
The energy functional E(h) is discretised using the same formalism, to give the following discrete energy:  . Input images before and after relative radiometric correction (Hessel et al., 2020) , normalised images and predicted heights from our shape-from-shading method for the first site (2.58 × 1.58 km). . Input images before and after relative radiometric correction (Hessel et al., 2020) , normalised images and predicted heights from our shape-from-shading method for the second site, more challenging because smaller (1.17 × 0.46 km).
this is a positive-definite quadratic form on the variable h. The unique minimum is attained at the point h solution of (N α Nα− λ∆)h = Nαf . Notice that, apart from the regularisation term depending on λ, this is just the normal equation of the first-order discrete problem. The proposed method consists in the solution of the regularised problem by this equation. The whole method can be implemented in about 10 lines of Matlab/Octave code.
So far, we have explained how to solve the linear shape-fromshading equation, which is just an approximation. To deal with the nonlinear problem we use the following strategy: first we solve the linear problem to obtain an initial approximation h 0 . Then we re-weight the data term f on the linear problem by the function 1 + (h 0 x ) 2 + (h 0 y ) 2 , which results in a new linear problem of the same form that has a unique solution h 1 . Thus we compute a sequence h n of approximate solutions that, if it converges, the limit must be the solution of the nonlinear problem. In practice, it seems that after a couple of iterations the solution does not change much.

EXPERIMENTS
For our experiments we use PlanetScope (Planet Team, 2017) imagery. The PlanetScope constellation is made of approximately 130 small satellites (form factor of 10 × 10 × 30 cm) imaging the entire Earth's landmass every day. The satellites carry 3-band or 4-band frame cameras and fly at 475 km on sunsynchronous orbits whose constant local solar time is between 9:30 and 11:30 am. The provided images are orthorectified with a GSD resolution of 3m.
Our dataset consists of a pair of temporal series of two sites over seven months, which amounts to about 100 images per site after removing cloudy days. We do not have access to the ground truth for the sites we studied. Still we were provided DEM reconstruction of the first site using stereo photogrammetry (de Franchis et al., 2014) obtained from SkySat imagery along with volumes (in cubic meters) computed from these DEMs. These reconstructions will serve as references for our validation.

Qualitative evaluation
We applied our method to the time series on the two sites and we observed that the stockpiles were correctly estimated in most cases. The images for the first site ( Figure 4) cover an area of 2.58 × 1.58 km while the images for the second site (Figure 5) only cover an area of 1.17 × 0.46 km. Thus the second site is much more challenging and we don't expect the same quality of results for both sites. Indeed our method struggles in reconstructing the shapes of the second site that presents much smaller stockpiles. By looking at the enhanced images we can confirm that it is actually quite hard to see the stockpiles in the images.

Analysis of time series
The ground is lower than the rails along which the cranes move. During the reconstruction we assumed the height of these rails to be zero. The first step to compute volumes is to add a global additive constant to restore the ground level to zero. The vertical scale of our SFS reconstruction is arbitrary but coherent in the whole time-serie. We add a multiplicative constant to this scale such that the maximum altitude of our stockpiles matches the one from the reference stereo reconstructions. The volumes (a) Colorbar Shape-from-Shading (b) Colorbar stereo Figure 6. Volume tracking in a time series. The plot represents the estimated volumes using the algorithm described in Section 3 on a time series of PlanetScope images spanning 7 months. The Shape-from-Shading heights are scaled by a constant so that the maximum altitude of our stockpiles matches the one from a reference stereo reconstruction (about 15m). The results for some highlighted dates are shown above the plot. The results of our method are consistent with the volumes obtained using stereo reconstruction and are obtained with a much higher frequency. The images shown below the plot allow the comparison on 3 dates of our result with 3D reconstructions. Note that the altitudes in the 3D reconstructions are expressed relative to the ellipsoid.
are then computed from the reconstructed shapes by integration on the ROI after clipping the first and last percentile to get rid of aberrant values. Then we multiply this result by 9, which is the size of a pixel for a GSD resolution of 3m.
In Figure 6 we represent the volumes estimated using the algorithm described in Section 3 on a time series of PlanetScope images spanning 7 months. The results for some highlighted dates are shown above the plot. We can see that the estimated volume is lower when less stockpiles are detected. For the stereo reconstructions we only have access to three months of data with about one acquisition per week. The images shown below the plot allow the comparison on 3 dates of the result of the 3D reconstruction with our method. Note that the PlanetScope and the SkySat images are often acquired at different hours or days. So the volumes obtained by stereo reconstruction are only an approximation of the ground truth. Taking this into account we see from Figure 6 that our estimated volumes are quite realistic.

CONCLUSION
Nowadays stockpile monitoring appears to be performed almost exclusively using UAV based techniques. In this paper we have shown that applying a simple Shape-from-Shading method on low resolution PlanetScope images can be sufficient to monitor the evolution of the volume of these stockpiles. The unquestionable loss in accuracy is well compensated in our opinion by the much higher revisit frequency and less human intervention. It is also easy to follow simultaneously several stocking sites all over the world. Thus this method appears very well suited to the problem of large scale stocks monitoring. In future work, more sophisticated Shape-from-Shading algorithms should be considered and the detection of occluding objects could be improved. Our method needs to be tested on more stocking sites. Access to more stereo reconstructions would also allow a quantitative evaluation of our results in addition to our qualitative evaluation.