Dense Local Azimuth-Elevation Map for the Integration of GIS Data and Camera Images

1. Introduction

We consider a situation where a pan-tilt-zoom (PTZ) camera fixed in position is used to monitor a landscape. In such a monitoring case, there is often an interest for associating to each camera pixel the geodetic coordinates (e.g. WGS84) [1,2] of the observed physical point. Performing this association is the georegistration of the image or the georeferencing [3] of the image pixels. With the geodetic coordinates of an observed point come along a plethora of possibly interesting geographical informations delivered by a Geographical Information System (GIS) [4].

Digital Elevation Models (DEM) provide information on the three-dimensional shape of the Earth’s surface in the form of elevation values (vertical datum) on a GIS raster (horizontal datum) [5]. Several governmental mapping agencies provide this type of topographical data. The Swiss Federal Office of Topography swisstopo e.g. provides the swissALTI3D [6] raster DEM. DEM are even available at the continent level for Europe as the EuroDEM dataset [7].

Regarding geometry, a camera performs a transformation from the 3D geometric space representing the physical world to the 2D space representing the image. This transformation has been widely studied in the fields of photogrammetry [8], computer vision [9,10], and computer graphics [11]. In most cases of modelling a real camera, this transformation is, in first approximation, a true perspective projection [9]. This is the often called pinhole camera model, referring to an ideal camera where all light rays pass through a pinhole. The transformation can be seen as being composed of three concatenated geometric transformations [9]. The coordinates of the source point in a 3D Cartesian world coordinate system are first transformed into its coordinates in a similar but camera-centric coordinate system, having the Z-axis aligned with the optical axis of the camera. This transformation depends on the camera extrinsic parameters: its pose, made of its position and its orientation. Then, occurs the projection to the image plane, which is perpendicular to the optical axis, giving rise to dimensionless 2D coordinates. Finally, these last coordinates are transformed into pixel coordinates by a 2D transformation depending on the camera intrinsic parameters: focal length, pixel size(s), and origin of the sensor coordinate system.

Except for wide-angle cameras [12], real cameras can be modelled with only a minor modification of the pinhole model [9]: a supplementary transformation is introduced between the second and third transformations mentioned above, in order to take into account the lens distortions. It consists in transforming the dimensionless coordinates resulting from the projection into distorted coordinates by adding a distortion term. It is common from at least the 1980’s [13] to decompose the distortions into a radial component (pure radial displacement) and a tangential component. One reason is that the radial component is often more important than the tangential one and this last one can be therefore often neglected. Furthermore, radial and tangential displacements, increasing with the distance from the image centre, are modelled by polynomial functions of the normalised coordinates.

Let us now assume that the model parameters of a camera monitoring a landscape are known and a DEM of this landscape is available. Then, georeferencing the image pixels is similar to rendering the surface spanned by the DEM with the camera from which we have the model. Indeed, in both cases a point or a small patch of the observed surface is assigned to each pixel. In the case of surface rendering, a colour is then attributed to the pixel according to the light emitted or re-emitted by the surface patch, while, in the case of georeferencing, the geodetic coordinates of the observed surface point or patch are attributed to the image pixel. Solutions for the rendering exist, in particular [14], and can be adapted for georeferencing.

Hence, in the case where a DEM of the landscape monitored by the camera is available, georeferencing the image pixels mainly boils down to finding the parameters of the camera model. This problem, known as geometric camera calibration, has been widely studied in the fields of photogrammetry [15,16] as well as computer vision [9,17,18] and a large number of methods have emerged in both communities. All have in common that the model parameters are estimated from a set of geometric primitives in the physical world and their corresponding geometric entities in the image. In probably the majority of the methods, the world geometric entities are points at the surface of physical objects and their corresponding entities are image pixels or, for more precision, points with sub-pixel coordinates. Prior to the computation of the parameter estimates itself, the first problem to be solved is the selection of calibration points, i.e. points in the scene and points in the image which correspond to each other.

There are two main approaches to solve the problem mentioned hereinbefore. The first one consists in introducing in the scene individual marks or patterns to materialise points, the corresponding points to which are easily detected in the image. Individual marks can be disks with a binary code allowing to identify each of them [19,20], or spheres [21,22]. A rigid planar chequerboard of known square size is used in the well known calibration method developed by Zhang [9,23,24], and implemented in the computer vision library OpenCV [25] as well as in the MATLAB^® Computer Vision Toolbox [26]. In this case, the selected points are the corners common to four of the squares.

The second approach is used in camera self-calibration [9,10,27,28,29], which involves multiple slightly different views of the same scene with the same camera. Points of interest are selected in the images and their correspondence from one image to the other is established based on the fact that they match, i.e. their neighbourhoods are similar. With the image points from the different views which correspond to each other, a virtual 3D point is constructed. It is the model of the physical point from which the image points originate. The correspondence of the image points with this virtual 3D point is in this way established automatically.

In the case we are interested in, i.e. a camera monitoring a landscape from which a DEM is at hand, none of the correspondence establishing methods mentioned above can be applied. For the first method, one cannot even imagine modifying the landscape, and accordingly the DEM, with marks detectable in the camera images. On the other hand, if it is easy with a PTZ camera, to generate multiple views of the same landscape, the generated virtual 3D points, mentioned above, will be referenced in a proper coordinate system and will therefore not be directly localisable in the DEM.

Selecting automatically calibration points in the image-DEM pair is the specific problem, this paper is concerned with. Moreover, the physical context is a camera having an oblique view on a landscape, i.e. far from an almost vertical top-down view. This article however does not go so far as to provide a complete solution to the problem. It proposes and describes a particular representation of the DEM, called local azimuth-elevation map (LAEM), the purpose of which is to ease the automatic selection of calibration points. A part from this interest, the LAEM could even be used in a second time for the image georegistration itself, in place of rendering the DEM with the camera model.

Section 2 is dedicated to related work. Then, in Section 3, the LAEM is presented as a particular 3D to 2D mapping. In the same section, a general method is described to compute an LAEM from a DEM in geodetic coordinates. Section 4 is concerned with implementation and experimental results. The implementation adapts the general method to the specificities of the swissALTI3D [6] DEM. Finally, Section 5 summarises the advantages of the particular DEM representation made by the LAEM and gives cues on how correspondences between an image and the LAEM could be found.

2. Related Work

To the best of my knowledge, the problem of selecting automatically points in the camera image and points in the DEM which correspond to each other, has been rarely addressed. However, since the motivation of this focused paper is image georegistration of oblique views, this section will present the related work in this broader view, with emphasis on the determination of the camera pose in a geographical coordinate system. Note that the latter is sometimes referred as the geolocalisation of the image, although the geolocalisation is often understood as concerning only the camera position without its orientation [30]. In the following, unless otherwise specified, geolocalisation will mean both position and orientation of the camera. The related work is organised according to the data used as reference for the georegistration, respectively the geolocalisation:

• orthophoto or georeferenced aerial images,
• georeferenced 3D point cloud,
• digital elevation model (DEM).

Instead of covering numerous publications, the overview gives priority to deep presentation of some examples in each category. Occasionally, personal comments are added.

Geo-Registration/Localisation with Orthophotos or Georeferenced Aerial Images.

Selecting corresponding points in the camera image and in an orthographic satellite image for calibrating the camera is performed semi-automatically in the software developed by Rameau et al. [31]. However, there are two limiting assumptions: the camera must have a rather top-down view on the scene and the calibration points have to lie in a more or less planar surface of the scene. The selection is semi-automatic, since the software user has first to rotate the satellite image until this one is almost aligned with the camera image. Then, a region of interest (ROI) has to be drawn, where the calibration points will be selected automatically.

Shan et al. [32] propose a solution for the georegistration of ground images with the help of georeferenced aerial images. They assume that the ground images have a GPS tag and hence an approximate geoposition. The method they propose is for improving the precision of the geolocation, adding also orientation to the position. Since, there are several ground images of the same scene, a self-calibration of the camera is made and, at the same time, a 3D model of the scene. With the self-calibration, the camera intrinsic parameters are known, but not its pose in a geographic coordinate system. The aerial views are not exactly vertical as is an orthophoto, but the difference between the angle of views of both type of images (terrestrial and aerial) is sufficiently large to make a direct matching between images unsuccessful. With the ground image, the 3D model, and the approximate geolocation of the ground image, they build a synthetic image from the aerial point of view. Then, by matching points between the synthetic image and the aerial image they determine the camera pose in the geographical coordinates system and hence the wanted more precise geolocalisation of the ground image. Comment. Shan et al. have shown impressive results with ground images of landmarks in Rome and aerial images of the city. These are images of man made shapes having a rather simple geometry. The method could be less successful with natural or rural landscapes.

Geo-Registration/Localisation with Georeferenced 3D Point Clouds.

Li et al. [33] propose a method for the geolocalisation of images with the help of georeferenced 3D point clouds. This is an extension to geolocalisation of previous work [34] or similar work by Sattler et al. [35], where the localisation was in a restricted area, typically a town. Let’s start with the localisation in a restricted area. What is called 3D point is in fact not a pure geometric primitive represented by coordinates in a 3D coordinate system, but features also descriptors taken from the images used to construct the geometric point. Hence, for establishing the correspondence between pixels in the image to be localised and some 3D points of the cloud, they use basically the same methods as for matching pixels between images. However, in creating constraints, the existing geometric information, helps for selecting the right corresponding 3D points in the vastness of the cloud. With the corresponding image pixels on one side and 3D points on the other side, the image can be localised by a camera calibration method. In the more recent work [33], the 3D point cloud has a worldwide extension, gathering 3D points from many places (cities) all around the world. A first difference with respect to the previous work is the significant larger size of the cloud, which makes the search of corresponding 2D-3D point pairs more difficult and necessitates a more performant solution. Another difference is that the localisation has to be done with respect to a global coordinate system, i.e. the image geolocalisation is now requested. They perform the georeferencing of the 3D points in using the geotag (latitude and longitude) of the images with which the 3D point has been constructed. As in the case of the localisation in a restricted area, the whole camera pose (3D position and orientation) is determined, but the geolocalisation seems to concern only the horizontal coordinates of the position, since they use precisely geotagged photos to test the precision of their geolocalisation.

DEM- Based Geo-Registration/Localisation.

Härer et al. [36] as well as Milosavljević et al. [37] propose solutions for georeferencing images with a DEM. However, for the step of selecting calibration points, they do it with the help of an orthophoto accompanying the DEM. Moreover, the selection is performed by a human operator. Comment. The manual selection of calibration points is time consuming and can even be difficult, in particular with natural landscapes. For example, ridge summits and passes selected in a camera image can be difficult to identify and localise precisely on an orthophoto or even on a topographical map. Proceeding in the opposite way, singular points detected in the orthophoto can be difficult to identify and localise in the camera image. This occurs when and why the view of the camera on the landscape is oblique, i.e. quite a part from the exactly vertical view of an orthophoto.

Portenier et al. [38] propose a solution for the georegistration of camera images with a DEM for the purposes of measuring snow cover. The camera position is assumed to be known a priori with sufficient precision, its orientation however not. Since they use a simple camera model, without lens distortions and with the optical image centre being in the middle of the digital image, only four parameters have to be estimated for calibrating the camera. Three are extrinsic: the orientation components, and one is intrinsic: the field of view (FOV), which is equivalent to the focal length presented in the introduction. The calibration process is automatic except for the first step, which consists in selecting at least one image with good sky-ground contrast, called the master image. In this image, the visible horizon is then detected and the polyline representing the horizon taken as reference. Then, the space of parameters is discretised within potential ranges of values. Next, for each discrete value of the parameter vector, the DEM is rendered with the camera model having these parameters, the visible horizon is detected and compared with the reference. Finally, the parameters selected for modelling the camera are those minimising the difference with the reference horizon. Portenier et al. noticed that the calibration could not be performed once and for all. Changes in weather conditions and/or human interactions, intentional or not, induce from time to time small changes in camera orientation, requiring a model update. Because it is a cumbersome process, the model is not updated by performing again the calibration. A matching is performed between an image of the moved camera and the master image. The camera rotation and hence the new orientation is then inferred from the pairs of corresponding image points.

Final Comments

The main difficulty in the geolocalisation with orthophotos and georeferenced aerial images is that the image to be geolocalised can differ significantly from the reference image(s). This is not only due to the difference of views, but also to the difference of capture times. In the time interval between the two captures (reference image and image to be geolocalised), drastic changes may have occurred in the scene (e.g. land coverage) or in its illumination (e.g. day/night). Because of these drastic changes, the images can be so different that the geolocalisation becomes impossible.

Geolocalisation with georeferenced 3D point clouds starts in fact with a large dataset of georeferenced photos taken by individuals, from which a georeferenced 3D point cloud is first generated. The 3D points serve to determine the camera pose of the image being geolocalised, but the correspondence between image pixels and some of the 3D points relies on the similitude of this image with the images that have been used to build the 3D points. Hence, the geolocalisation with georeferenced 3D point clouds is subject to the same problem of difference in capture time between images as mentioned before. However in this case, the problem can be mitigated in having at hand many images of the same scene, taken under different lighting and weather conditions. The fact that the 3D point clouds stem from photos taken by individuals during their leisure time has the effect of restricting the use of this geolocalisation method to places which are popular and very much visited, e.g. landmarks in cities.

DEM- based geolocalisation can be associated to geolocalisation with georeferenced 3D point clouds, since a DEM is indeed a georeferenced 3D point cloud. However, the physical reality modelled by the point cloud differs. As seen in the previous paragraph, in the approach named “with 3D point clouds”, the 3D geometry is the one of very much visited built spaces, like landmarks in cities. DEMs in contrary are modelling the ground or, scientifically speaking [5], the interface between the atmosphere and the lithosphere, of large regions. The approaches are therefore not at all rivals but complementaries.

Another difference between both approaches is that DEMs are not built with images from multiple views and therefore do not hold image information which can be used to find the correspondences between the image to be geolocalised and the DEM. It is true that there often exists an orthophoto which can be associated point by point with the DEM. However, because of the problematic large difference in views mentioned in the paragraph dedicated to the geolocalisation with orthophotos, the similarity between the orthophoto and the image to be geolocalised is poor. The operator selecting the calibration points in Milosavljević et al. [37] undoubtedly uses his understanding of the image to be geolocalised and of the orthophoto to make the selection. New paths have to be open for an automatic selection of corresponding points between the image to be geolocalised/georegistered and a DEM. Such an attempt has been made by Portenier et al. [38]. This paper proposes a new one.

3. Local Azimuth-Elevation Map at Camera Position

3.1. Scenario

As explained in the introduction, we consider the case of a PTZ camera fix in position with an oblique view on a natural landscape from which we have a DEM. The typical camera set-up, which is the one used for the experiments (Section 4), is represented in Figure 11. The camera pan axis is more or less vertical and the pan rotation allows a 360° view around the camera position. The proposed method is not exclusive to PTZ cameras. But it is more useful in such a case, where the camera has changing orientation and changing field of view, and it is especially well suited for the camera set-up shown in Figure 1.

Since the camera has a fix position, determining its geodetic coordinates is a one-time procedure, which can be performed at the camera installation or even prior to it. It is is an easy task, especially if the observed landscape is relatively distant from the camera, since in this case, the required positioning accuracy is low. The camera geodetic coordinates can be determined e.g. in using a hand-held global positioning system (GPS). They can also be read on a digital map e.g. on Google Maps [40].

Knowing a priori the camera position and assuming that the DEM is available in the same geodetic coordinate system as the camera position , we will see hereafter that we can represent the DEM in a local coordinate system, where the origin is at the camera position.

3.2. The ENU and AER Local Coordinate Systems

You can refer to Figure 2 for a graphical support to the following text.

Let us consider a point $P_0$, the camera position in our case, with geodetic coordinates [1,2,41]: geodetic latitude ${\varphi }_0$, longitude ${\lambda }_0$, and geodetic or ellipsoidal height $h_0$. Let $S_0$ be the orthogonal projection of $P_0$ on the ellipsoid modelling the Earth. $h_0$ is hence the Euclidian norm of the vector $\overrightarrow{S_0{P}_0}$. Let ${\pi }_0$ be the plane tangential to the Earth ellipsoid at $S_0$ and $\pi$ be the plane parallel to ${\pi }_0$ and containing $P_0$.

The local East-North-Up (ENU)[1,41] or Local Geodetic Horizon (LZH) [2] coordinate system with origin $P_0$ is a Cartesian coordinate system having the three following to each other orthogonal axis:

• Up is the direction of the vector $\overrightarrow{S_0{P}_0}$, denoted by $z_U$ in Figure 2.
• North is the direction at $S_0$ of the meridian, i.e. the ellipse through the geodetic North. It is taken positive towards the geodetic North . In Figure 2, it is represented in the plane $\pi$ parallel to ${\pi }_0$, and denoted by $y_N$.
• East is the direction orthogonal to the two other ones and so that ENU is a right-handed system. It is denoted by $x_E$ in Figure 2.

Since ENU is a Cartesian coordinate system, it can be used to represent 3D points for camera calibration.

The formulas to transform the coordinates of a point in the geodetic coordinate system to its coordinates in the local ENU coordinate system use an intermediate representation in the Earth-centred, Earth-fixed (ECEF) coordinate system. They can be found e.g. in  [1,2]. They are implemented in the MATLAB^® Mapping Toolbox by the function geodetic2enu [42] and in the the PyMap3D Python module by the function pymap3d.enu.geodetic2enu [43].

The coordinate transformation is invertible, although there is no explicit form for the geodetic coordinates in function of the ECEF coordinates [1]. The transformation from the ENU coordinates into the geodetic coordinates is implemented in the MATLAB^® Mapping Toolbox by the function enu2geodetic [44] and in the the PyMap3D Python module by the function pymap3d.enu.enu2geodetic [45].

The Azimuth-Elevation-Range (AER) coordinate system is the spherical companion of the ENU system  [41]. Azimuth and elevation angles are a widely used coordinate system in astronomy to represent the direction from an observer position [46]. The elevation angle is also called sometimes altitude angle. There exist also some variations in the way these angles are measured. In the AER coordinate system, the coordinates of a point P also include the range (or slant range) i.e. the distance from $P_0$ to P. In the following, we use for the AER coordinate system the same conventions as in [41]. Let $P_{\pi }$ be the orthogonal projection of the point P on the plane $\pi$, then:

• the azimuth $\phi$ is the clockwise angle from the axis $y_N$ to the direction $\overrightarrow{P_0{P}_{\pi }}$;
• the elevation (or altitude) $\theta$ is the angle between the plane $\pi$ and the direction $\overrightarrow{P_{0}P}$, positive up, negative down;
• (slant) range r is the Euclidean distance between the points $P_0$ and P.

As in  [41], the angles will be expressed hereafter in degrees, in the range $[0^{\circ },{360}^{\circ }[$ for the azimuth $\phi$ and in the range $[-{90}^{\circ },{90}^{\circ }]$ for the elevation $\theta$.

The formulas for computing the AER coordinates from the ENU coordinates are:

$$\begin{array}{cc}\hfill \phi & =\left\{\begin{array}{cc}atan2(x_E,y_N)\hfill & \mathrm{if}{x}_E\ge 0,\hfill \\ atan2(x_E,y_N)+{360}^{\circ }\hfill & \mathrm{if}{x}_E<0,\hfill \end{array}\right.\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill \theta & =atan2(z_U,\sqrt{x_E^2+y_N^2}),\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill r& =\sqrt{x_E^2+y_N^2+z_U^2}.\hfill \end{array}$$

(1c)

Note that the azimuth $\phi$ is undefined when both $x_E$ and $y_N$ are null. The elevation $\theta$ is undefined when all three ENU coordinates are null, i.e. for the coordinate system origin $P_0$.

The function geodetic2aer [47] from the MATLAB^® Mapping Toolbox and the function pymap3d.aer.geodetic2aer [48] from the PyMap3D Python module allow to perform directly the whole transformation from the geodetic coordinates to the AER coordinates.

The formulas for computing the ENU coordinates from the AER coordinates are:

$$\begin{array}{cc}\hfill x_E& =rcos\theta sin\phi ,\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill y_N& =rcos\theta cos\phi ,\hfill \end{array}$$

(2b)

$$\begin{array}{cc}\hfill z_U& =rsin\theta .\hfill \end{array}$$

(2c)

The function aer2geodetic2 [49] from the MATLAB^® Mapping Toolbox and the function pymap3d.aer.aer2geodetic [50] from the PyMap3D Python module allow to perform directly the whole transformation from the AER coordinates to the geodetic coordinates.

3.3. Local Azimuth-Elevation Map (LAEM)

3.3.1. Mapping to the Azimuth-Elevation Space

The local azimuth-elevation map (LAEM) is defined as the result of a mapping from the 3D geometric space ${\mathbb{R}}^3$ to the 2D space ${\mathbb{R}}^2$ where the azimuth angle $\phi$ and the elevation angle $\theta$ of a source point are represented in a Cartesian coordinate system. Figure 3 represents this mapping graphically for a point $P_n$. Since the azimuth $\phi$ is undefined if both $x_E$ and $y_N$ are null, the points belonging to the $z_U$-axis do not have an image. Furthermore, the image of ${\mathbb{R}}^3$ is a limited domain of ${\mathbb{R}}^2$. With the conventions adopted in Section 3.2 for the azimuth and elevation angles, the image domain is given by the following ranges on $\phi$ and $\theta$ values: $(\phi \in [0^{\circ },{360}^{\circ }\left[\right)\wedge (\theta \in ]-{90}^{\circ },{90}^{\circ }\left[\right)$. Note that the domain does not include the elevation values $-{90}^{\circ }$ and ${90}^{\circ }$, because they correspond to undefined azimuth values (points on the $z_U$-axis).

Even if it is due to a well known property of the inverse trigonometric functions, let us point out the following feature of the image domain. Two points close to each other in the 3D space, one with a positive $x_E$ value and one with a negative $x_E$ value, have distant images. The first image point has a $\phi$ value slightly larger than $0^{\circ }$, while the second one has a $\phi$ value slightly smaller than ${360}^{\circ }$.

The hereinbefore introduced mapping is similar to the perspective projection realised by a camera of position $P_0$, in the sense that all points of a half-line from $P_0$ are mapped to a single point. However, there are at least two differences:

• the formula to derive the image coordinates $(\phi ,\theta )$ from the AER coordinates of the 3D point is different from the one for deriving the normalised camera image coordinates $(x,y)$ from the camera centric Cartesian coordinates $(X,Y,Z)$,
• a line segment, the endpoints of which map to different points, does not map into a line segment.

3.3.2. LAEM of a Surface

The LAEM of a surface in the 3D space is the mapping of the points belonging to the surface, however with an exclusion condition: from the points belonging to the same half-line from $P_0$, only the nearest from $P_0$, i.e. the "visible" one, is mapped. With this condition, the mapping becomes injective. As a consequence, for each pair of values $(\phi ,\theta )$ there exists maximum one surface point mapped to this image point. In Figure 4, the points $P_n$ and $P_n^{\prime }$ have the same azimuth ${\phi }_n$ and elevation ${\theta }_n$ values. However, only the point $P_n$ with smaller range value $r_n$ is mapped to the image point $({\phi }_n,{\theta }_n)$. As a consequence, the unique surface point $P_n$ can be associated to the values ${\phi }_n$ and ${\theta }_n$.

3.3.3. Discrete LAEM

For the purpose explained in the introduction, the LAEM has to be discrete, i.e. be a countable set of values. Moreover, the set has to be finite. Among the various possibilities of discretising the LAEM, let us propose this one. The choice is, at least partly, driven by the goal of representing a discrete LAEM as a digital image for an easy storage and processing.

We use a constant sampling interval $\Delta \phi$ in azimuth as well as a constant sampling interval $\Delta \theta$ in elevation. For convenience, we require that $\Delta \phi$, respectively $\Delta \theta$, divides ${90}^{\circ }$, formally:

$$\begin{array}{c}\hfill {\displaystyle \frac{{90}^{\circ }}{\Delta \phi }}=⌊{\displaystyle \frac{{90}^{\circ }}{\Delta \phi }}⌋,\end{array}$$

(3a)

$$\begin{array}{c}\hfill {\displaystyle \frac{{90}^{\circ }}{\Delta \theta }}=⌊{\displaystyle \frac{{90}^{\circ }}{\Delta \theta }}⌋.\end{array}$$

(3b)

In practice (see Section 4), this constraint is not restrictive, since we use $\Delta \phi$ and $\Delta \theta$ values as small as $0.{01}^{\circ }$.

As it is usual to index the rows of a digital image with i, starting at the top, and the columns with j, starting at left, we index the discrete azimuth values ${\phi }_j$ with j, starting with the smallest, and the discrete elevation values ${\theta }_j$ with j, starting with the largest. This is represented in Figure 5, where x and y, the commonly used image coordinates [51,52], are also represented. Most of the programming languages start table indexes with 0, this is the assumption for i and j. MATLAB^® in contrary starts table indexes with 1. We use the identifiers $i^{\prime }$ and $j^{\prime }$ for this case. The relationships between both types of indexes are of course simply given by:

$$\begin{array}{ccc}\hfill i^{\prime }& =i+1,\hfill & \hfill i=i^{\prime }-1,\end{array}$$

(4a)

$$\begin{array}{ccc}\hfill j^{\prime }& =j+1,\hfill & \hfill j=j^{\prime }-1.\end{array}$$

(4b)

In the envisioned camera to landscape configurations, only a limited portion of the elevation range $]-{90}^{\circ },{90}^{\circ }[$ is significant. Let ${\theta }_{\min }$ and ${\theta }_{\max }$ be the minimal, respectively the maximal, elevation value of interest. For convenience, we require ${\theta }_{\min }$ and ${\theta }_{\max }$ be both integer multiples of $\Delta \theta$, formally:

$$\begin{array}{cc}\hfill {\theta }_{\min }& =k_1·\Delta \theta ,k_1\in \mathbb{Z},\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill {\theta }_{\max }& =k_2·\Delta \theta ,k_2\in \mathbb{Z}.\hfill \end{array}$$

(5b)

With these conditions, the number $N_{\theta }$ of discrete elevation values ${\theta }_i$ is

$$N_{\theta }={\displaystyle \frac{{\theta }_{\max }-{\theta }_{\min }}{\Delta \theta }}+1$$

(6)

and the discrete elevation values ${\theta }_i$ are given by the equation

$${\theta }_i={\theta }_{\max }-i·\Delta \theta .$$

(7)

Concerning the azimuth, there may be cases where the full range of values $[0^{\circ },{360}^{\circ }[$ is of interest. However, there also may be cases where the interest is only for a segment of the possible range. It will of course occur that this segment of interest will be on both sides of the azimuth $0^{\circ }$. In the case where the full range of azimuth values is of interest, one could decide to handle easily the angle wrapping, i.e. the transition from $359.{\overline{9}}^{\circ }$ to $0^{\circ }$, in expanding the LAEM with an overlap corresponding to the horizontal field of view of the camera. In order to handle all these cases in a simple way, we propose to use unwrapped angle values, i.e. angles values always increasing. Let ${\phi }^{\prime }$ denotes the unwrapped azimuth and ${\phi }_{\min }^{\prime }$ and ${\phi }_{\max }^{\prime }$ the minimal, respectively, the maximal, unwrapped angle value. For the full range with overlap mentioned above, one would set for example ${\phi }_{\min }^{\prime }=-{\displaystyle \frac{\alpha }{2}}$ and ${\phi }_{\max }^{\prime }={360}^{\circ }+{\displaystyle \frac{\alpha }{2}}$, where $\alpha$ is the horizontal field of view of the camera . For convenience, we require ${\phi }_{\min }^{\prime }$ and ${\phi }_{\max }^{\prime }$ be both integer multiples of $\Delta \phi$, formally:

$$\begin{array}{cc}\hfill {\phi }_{\min }^{\prime }& =k_1·\Delta \phi ,k_1\in \mathbb{Z},\hfill \end{array}$$

(8a)

$$\begin{array}{cc}\hfill {\phi }_{\max }^{\prime }& =k_2·\Delta \phi ,k_2\in \mathbb{Z}.\hfill \end{array}$$

(8b)

With these conditions, the number $N_{{\phi }^{\prime }}$ of discrete unwrapped azimuth values ${\phi }_j^{\prime }$ is:

$$N_{{\phi }^{\prime }}={\displaystyle \frac{{\phi }_{\max }^{\prime }-{\phi }_{\min }^{\prime }}{\Delta \phi }}+1$$

(9)

and the discrete unwrapped azimuth values ${\phi }_j^{\prime }$ are given by the equation

$${\phi }_j^{\prime }={\phi }_{\min }^{\prime }+j·\Delta \phi .$$

(10)

The discrete wrapped azimuth angle ${\phi }_j$ corresponding to ${\phi }_j^{\prime }$ is then simply given by

$${\phi }_j=\mathrm{mod}({\phi }^{\prime },{360}^{\circ }),$$

(11)

where the mod function is defined by

$$\mathrm{mod}(a,b)=a-b·⌊{\displaystyle \frac{a}{b}}⌋.$$

(12)

3.4. Computing a Discrete LAEM of a Gridded DEM

3.4.1. Methods from Surface Rendering

Section 3.3.1 and Section 3.3.2 show up that computing a discrete LAEM of a surface in 3D space is close to rendering this surface with a camera model, a basic operation in computer graphics [11], and a DEM is basically a surface. It is therefore worth considering if methods used for surface rendering could be adapted to the computation of a discrete LAEM.

Surface rendering starts by establishing the correspondence between image points, i.e. pixels, and points belonging to the surface. For determining this correspondence, there are basically two opposite approaches. The first one proceeds from the scene to the image. The surface is modelled by polygonal patches, either as an exact model or as an approximation obtained by tessellation. Then, the image of the polygon vertices is computed by the true perspective projection modelling the camera. Finally, since linear segments are imaged into linear segments, the polygon edges are drawn as well as the polygon interior areas. This approach is less precise than the second one, because the correspondence is not established individually for each pixel within a patch. However, it is less computation intensive. It is therefore preferred, when the computation time is critical. Whatever the advantages of this approach, it cannot be adapted to LAEM computation, since in this case, as mentioned in Section 3.3.1, line segments are not mapped into line segments.

The second approach proceeds from the image to the scene. This method is known as ray casting since [53]. A ray is traced from the camera centre through each image pixel and the corresponding surface point is determined as the intersection of the ray with the surface. In the case of an LAEM, points of the 3D space of constant elevation and azimuth build also rays from a fixed point, here the coordinate system origin. The only difference with respect to ray casting lies in the distribution of rays in the space. In the case of surface rendering, rays are regularly spaced according to the angle tangent, while in the case of LAEM, they are regularly spaced with respect to the angle value itself. The ray casting method can therefore be easily adapted to LAEM computation. In this paper, we are merely mentioning this possibility, emphasising the method presented in the next section.

3.4.2. Proposed Method: Principle

The method, this paper promotes for computing a discrete LAEM, is different from the surface rendering methods discussed in the previous section. It proceeds mainly from the scene (DEM) to the image (LAEM), but also partly from the LAEM to the DEM. Moreover, it is not linked to rendering the surface represented by the DEM, since it alternatively makes use of the volume of which the DEM is the upper boundary.

Let $P_{m,n}$ be a particular element of the DEM with position indices m and n and let us take one of the 4-neighbourhood elements, e.g. $P_{m-1,n}$. The elementary surface patch used for the LAEM computation is a delimited cross-section of the volume associated with the DEM, as represented in Figure 6. The delimiting edges are: the line segment $[P_{m-1,n},P_{m,n}]$, the two segments of vertical lines from $P_{m-1,n}$, respectively $P_{m,n}$, until the minimal elevation value ${\theta }_{\min }$ is reached, and the line segment joining these two last points.

The motivation for taking vertical lines is that a DEM provides elevation values (vertical datum) on a GIS raster (horizontal datum) [5]. At the origin of the local ENU coordinate system, the DEM vertical direction is the one of the Up-axis. However, the further a DEM point is from the ENU origin, the larger is the deviation of the vertical direction at that point from the Up-axis of the ENU coordinate system. For a DEM where the elevation is the geodetic height, as we have assumed, the vertical direction is the one of the normal to the underlying ellipsoid. Determining the exact vertical direction at any point of a DEM is however cumbersome and there is generally no need to be very precise on the vertical direction. To have a reference value for the needed precision, let us take a DEM grid resolution of 2 meters, like the DEM used in the experimentation, and the greatest purely vertical drop on the Earth [54], which is 1200 m. The angle formed with respect to the vertical by two consecutive DEM points, one at the top and one at the bottom of the drop, is $arctan(2/1200)$, equal to $0.{0955}^{\circ }$. This, rounded to $0.1^{\circ }$, is the sufficient precision in determining the vertical direction. Given this tolerance, we propose to determine the vertical direction with one among the two following degrees of approximation of the DEM horizontal datum: planar and spheric (see Section 3.4.3).

As explained in the introductory paragraph of the section, the method also involves a procedure from the LAEM to the DEM. It is the following. Let ${\phi }_{j1}^{\prime }$ and ${\phi }_{j2}^{\prime }$ be the discrete azimuth values corresponding to $P_{m-1,n}$, respectively $P_{m,n}$, as represented in Figure 6. If the indices $j_1$ and $j_2$ are apart from more than one, i.e. formally: $|j_2-j_1|>1$, then for any in-between value j, i.e. $j_1<j<j_2$, a point $P_j$ is determined on the line segment $[P_{m-1,n},P_{m,n}]$ as the intersection of the plane defined by all the points of same azimuth ${\phi }_j^{\prime }$ and the line passing through the points $P_{m-1,n}$ and $P_{m,n}$ (see Figure 7).

The ENU coordinates of the point $P_j$ can be easily computed from ${\phi }_j^{\prime }$ and the ENU coordinates of $P_{m-1,n}$ and $P_{m,n}$ by solving a system of linear equations. Let $(x_1,y_1,z_1)$ be the coordinates of $P_{m-1,n}$ in the ENU coordinate system and $(x_2,y_2,z_2)$ be those of $P_{m,n}$. Then, the coordinates $(x_j,y_j,z_j)$ of $P_j$ in the ENU coordinate system can be found by solving the system of equations:

$$\begin{array}{cc}\hfill -cos{\phi }_j^{\prime }·x_j+sin{\phi }_j^{\prime }·y_j& =0,\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill x_j& =\lambda (x_2-x_1)+x_1,\hfill \end{array}$$

(13b)

$$\begin{array}{cc}\hfill y_j& =\lambda (y_2-y_1)+y_1,\hfill \end{array}$$

(13c)

$$\begin{array}{cc}\hfill z_j& =\lambda (z_2-z_1)+z_1.\hfill \end{array}$$

(13d)

(13a) is the equation of the plane determined by ${\phi }_j^{\prime }$, while equations (13b) to (13d) are the parametric equations of the line through the points $P_{m-1,n}$ and $P_{m,n}$. Equation (13a) can be replaced by the following explicit equation for $\lambda$:

$$\begin{array}{cc}\hfill \lambda & ={\displaystyle \frac{cos{\phi }_j^{\prime }·x_1+-sin{\phi }_j^{\prime }·y_1}{-cos{\phi }_j^{\prime }·(x_2-x1)+sin{\phi }_j^{\prime }·(y2-y1)}}.\hfill \end{array}$$

(14a)

Once the point $P_j$ has been determined, the patch is "filled" in following down the vertical direction from this point.

3.4.3. Computing DEM vertical lines

As explained in Section 3.4.2, for computing vertical lines from a DEM point, we propose to use approximations of the DEM horizontal datum. We consider two possible approximations: planar and spheric.

Planar Horizontal Datum 

If the DEM horizontal datum is the local horizontal plane, as represented in Figure 3, the DEM vertical direction is the normal to this plane. Hence, all points on a vertical line passing through a DEM point have the same coordinates $x_E$ and $y_N$ as this point, as well as the same azimuth $\phi$ and horizontal distance $\varrho =\sqrt{x_E^2+y_N^2}$.

Spheric horizontal datum. 

At $P_0$, the origin of the local coordinate system, the curvature of the ellipsoid underlying the DEM depends on the direction considered. According to [1], the minimal radius of curvature is in the meridian, i.e. in the $y_N$-axis direction, while the maximal is in the perpendicular direction, i.e. the one of the $x_E$-axis. This last curvature is named curvature in the prime vertical. The values of these two extreme local radii of curvature depend on the geodetic latitude ${\varphi }_0$ of $P_0$. The formulae for the radius of curvature in the meridian $\rho \left({\varphi }_0\right)$ and for the radius of curvature in the prime vertical $\nu \left({\varphi }_0\right)$ are:

$$\begin{array}{cc}\hfill \rho \left({\varphi }_0\right)& ={\displaystyle \frac{a(1-e^2)}{{(1-e^2{sin}^2{\varphi }_0)}^{3/2}}},\hfill \end{array}$$

(15a)

$$\begin{array}{cc}\hfill \nu \left({\varphi }_0\right)& ={\displaystyle \frac{a}{{(1-e^2{sin}^2{\varphi }_0)}^{1/2}}},\hfill \end{array}$$

(15b)

where a is the length of the ellipsoid major axis, 6378137.0 m for GRS80, and $e^2$ is the ellipsoid eccentricity squared, 0.006694380022901 for GRS80.

As an approximation of the DEM horizontal datum around $P_0$ by a sphere, we propose to use a radius R between the radius of curvature in the meridian $\rho \left({\varphi }_0\right)$ (minimal value) and the radius of curvature in the prime vertical $\nu \left({\varphi }_0\right)$ (maximal value), of course augmented by the geodetic height $h_0$ of $P_0$. Since the geometric mean of both curvatures has a simpler formula than the arithmetic mean, we give preference to the first option. The formula for computing R from the geodetic latitude ${\varphi }_0$ and height $h_0$ is:

$$\begin{array}{cc}\hfill R& =\sqrt{\rho \left({\varphi }_0\right)\nu \left({\varphi }_0\right)}+h_0,\hfill \end{array}$$

(16a)

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{a\sqrt{1-e^2}}{1-e^2{sin}^2{\varphi }_0}}+h_0.\hfill \end{array}$$

(16b)

Given a DEM point $P_n$, the vertical line from this point is normal to the sphere and hence crosses the $z_U$-axis at the sphere centre C. Let us consider a plane containing $P_n$ and parallel to the horizontal plane through $P_0$. We can draw a rectangle in this plane as represented in Figure 8(a). The lengths of its edges are the coordinates $x_E$ and $y_N$ of $P_n$. For any point on the vertical line from $P_n$, we can do the same: take a parallel plane and draw the rectangle. The edge lengths of each rectangle are the coordinates $x_E$ and $y_N$ of the corresponding point on the vertical line. All the rectangles are similar to each other. Consequently, the coordinates $x_E$ and $y_N$ of the points on the vertical line have a constant ratio $x_E/y_N$ or $y_N/x_E$, and therefore a constant azimuth, equal to the one one $P_n$, that is ${\phi }_n$.

Figure 8(a) represents the vertical line from $P_n$ in the plane $C{P}_0{P}_n$. Considering the triangle $C{P}_n{P}_U$, it is obvious that the angle $\alpha$ between the vertical line from $P_n$ and the $z_U$-axis depends on the elevation angle ${\theta }_n$ and the range $r_n$ of $P_n$ according to the following formula:

$$tan\alpha ={\displaystyle \frac{r_{n}cos{\theta }_n}{R+r_{n}sin{\theta }_n}}.$$

(17)

The elevation and the range of any point on the vertical line from $P_n$ satisfy the same equation, in particular the points corresponding to the discrete elevation values ${\theta }_i$. The implicit equation for ${\theta }_i$ and the corresponding range $r_i$ can be rewritten into:

$$r_i(cos{\theta }_i-tan\alpha sin{\theta }_i)=Rtan\alpha ,$$

(18)

which can be converted into the following explicit equation for $r_i$:

$$r_i={\displaystyle \frac{Rtan\alpha }{cos{\theta }_i-tan\alpha sin{\theta }_i}}.$$

(19)

Maximal Errors in Vertical Direction 

We end this section by proposing an estimated upper bound for the errors in vertical direction committed with each of the two horizontal datum approximations presented hereinbefore.

Let us consider a point $P_n$ in the local horizontal plane distant from $P_0$ by d. This is a particular case of Figure 8(b). The vertical line from $P_n$ is for a spheric model or radius R. The equation (17) for determining the angle $\alpha$ between the vertical line from $P_n$ and the $z_U$-axis takes the simpler form: $tan\alpha =d/R$.

In the case of the planar approximation, the vertical direction at $P_n$ given by the model is the one of the $z_U$-axis. If we assume that the right vertical direction is the one given by the sphere of radius R, the error is measured by the angle $\alpha$. For a given distance d, ${\alpha }_R$ is maximal for R being minimal. The minimal value for R is the radius of curvature in the meridian at the equator, i.e. $\rho \left(0^{\circ }\right)$. We hence consider ${\alpha }_{\rho \left(0^{\circ }\right)}$ as an upper bound of the error in vertical direction made in using the planar approximation.

Around $P_0$, the ellipsoid modelling the Earth lies between a sphere of radius $\rho \left({\varphi }_0\right)$ and a sphere of radius $\nu \left({\varphi }_0\right)$. The angle error in estimating the vertical direction with a sphere of mean radius, as we do with the spheric approximation, should be smaller than the difference between the vertical directions determined by both spheres. This difference is maximal at the equator. We consider therefore $({\alpha }_{\rho \left(0^{\circ }\right)}-{\alpha }_{\nu \left(0^{\circ }\right)})$ as an upper bound of the error in vertical direction made in using the spheric approximation.

Table 1 gives the numerical values of ${\alpha }_{\rho \left(0^{\circ }\right)}$, respectively $({\alpha }_{\rho \left(0^{\circ }\right)}-{\alpha }_{\nu \left(0^{\circ }\right)})$ at different distances d of $P_0$. The used reference ellipsoid is GRS80. Given the tolerance of $0.1^{\circ }$ for a grid resolution of 2 meters computed in Section 3.4.2, it follows from the numerical values of Table 1 that the planar approximation can be used up to a distance of about 10km and that the spheric approximation has no practical limits for this tolerance.

3.4.4. Detailed Method Description

We use a scalar two-dimensional array to represent the LAEM. Hereafter, we will denote an array element by $d_{i,j}$ and the whole array by $\mathbf{D}$. In the case of the planar approximation, d is the horizontal distance $\varrho =\sqrt{x_E^2+y_N^2}$, while in the case of the spheric approximation, d is the slant range r. In any of these cases, the triplet $({\phi }_j,{\theta }_i,d_{i,j})$ represents a point of the 3D space, the coordinates of which can be computed in any coordinate system: AER, ENU, geodetic,…

The pseudo-code of listing 1 explains in detail how the LAEM $\mathbf{D}$ is computed from a gridded DEM, the elements of which are denoted by $P_{m,n}$. We only add information where we estimate necessary, referring to the code by the line number.

Listing 1. LAEM computation
1	% compute coordinates of DEM points $P_{m,n}$
2	for each $P_{m,n}$:
3	compute $j_{m,n},i_{m,n},d_{m,n}$
4	% add coordinates of interpolation points $P_j$
5	for each $j_{m,n}$:
6	for each $j\mid j_{m,n}<j<j_{m+1,n}$ and each $j\mid j_{m,n}<j<j_{m,n+1}$:
7	compute $i_j,d_j$ making use of (14a), (13b), (13c), and (13d)
8	% fill the LAEM $\mathbf{D}$
9	$d_{i,j}=+\infty ,\forall i,j$
10	for each $(j_P,i_P,d_P)\in \left\{(j_{m,n},i_{m,n},d_{m,n})\right\}\cup \left\{(j,i_j,d_j)\right\}$:
11	if $d_P<d_{i_P,j_P}$:
12	for each $i\in \{i_P,i_P+1,\dots ,N_{\theta }-1\}$:
13	% spheric approximation	%planar approximation
14	compute $r_i$ with equation (19)	nop
15	$d_{i,j_P}=r_i$	$d_{i,j_P}=d_P$

line 3: 

The values $j_{m,n}$, $i_{m,n}$, and $d_{m,n}$ result directly from the representation of $P_{m,n}$ in the AER coordinate system. $d_{m,n}$ represents either the horizontal distance (planar approximation) or slant range (spheric approximation) of $P_{m,n}$. $j_{m,n}$ and $i_{m,n}$ determine the position of $P_{m,n}$ in the LAEM. They are derived from the azimuth $\phi$ and the elevation $\theta$ of $P_{m,n}$ according to the choices that have been made in Section 3.3.3 for the discretisation of the $(\phi ,\theta )$ space. The equations are:

$$\begin{array}{cccc}\hfill j& =⌊{\displaystyle \frac{{\phi }^{\prime }-{\phi }_{\min }^{\prime }}{\Delta \phi }}+{\displaystyle \frac{1}{2}}⌋,\hfill & \hfill & \mathrm{if}{\phi }_{\min }^{\prime }-{\displaystyle \frac{\Delta \phi }{2}}<{\phi }^{\prime }\le {\phi }_{\max }^{\prime }+{\displaystyle \frac{\Delta \phi }{2}},\hfill \end{array}$$

(20a)

$$\begin{array}{cccc}\hfill i& =⌊{\displaystyle \frac{{\theta }_{\max }-\theta }{\Delta \theta }}+{\displaystyle \frac{1}{2}}⌋,\hfill & \hfill & \mathrm{if}{\theta }_{\min }-{\displaystyle \frac{\Delta \theta }{2}}<\theta \le {\theta }_{\max }+{\displaystyle \frac{\Delta \theta }{2}},\hfill \end{array}$$

(20b)

where ${\phi }^{\prime }$ is the unwrapped azimuth angle, related to the azimuth angle $\phi$ by:

$${\phi }^{\prime }=\phi +k·{360}^{\circ },k\in \{-1,0,1\}.$$

(21)

line 7: 

The interpolation has been explained in Section 3.4.2 and illustrated in Figure 7. j determines the discrete unwrapped azimuth angle ${\phi }_j^{\prime }$. The explicit equations (14a), (13b), (13c), and (13d) allow to compute from ${\phi }_j^{\prime }$ the ENU coordinates of the interpolation point $P_j$ between $P_{m,n}$ and $P_{m+1,n}$, respectively between $P_{m,n}$ and $P_{m,n+1}$. From the ENU coordinates of the interpolation point, the corresponding $\theta$ and $d_j$ values are derived with the ENU to AER coordinate transform. Finally, $i_j$, determining together with j the position of $P_j$ in the LAEM, is computed from $\theta$ by equation (20b). This new coordinate triplet $(j,i_j,d_j)$ is added to the set (or list) of coordinates that will be browsed to fill the LAEM.

lines 

9 and 11: According to Section 3.3.2, from the points having same LAEM position, we keep only the "visible" one, i.e. the one with smallest d value. The LAEM is hence initialised at line 9 with a specific value, larger than any possible d value, here the choice $+\infty$. Furthermore at line 11, the position $(i_P,j_P)$ of the currently evaluated point is only updated, if the value $d_P$ is smaller than the value stored in $\mathbf{D}$ at this position.

lines 

12 to 15: $i_P$ together with $j_P$ determine the position of the currently evaluated DEM point. The following i values, each together with $j_P$, determine the LAEM positions of 3D points located on the vertical line down from the DEM point. According to Section 3.4.3, they have the same azimuth and hence the same j value. Moreover, in the case of the planar approximation, they have the same horizontal distance $\varrho$ value as the as the DEM point, i.e. $d_P$. In the case of the spheric approximation, the point slant range $r_i$ is computed from ${\theta }_i$ by equation (19). The LAEM is updated with these respective values: $d_P$ and $r_i$. Note that it is not necessary to test if the d value of the points on the vertical line is smaller than the value existing in the LAEM at that position, because the truth follows directly from the truth of the condition for the DEM point.

3.4.5. DEM Downsampling as Preprocessing

The discrete LAEM as defined in Section 3.3.3 has a constant azimuthal resolution $\Delta \phi$. The corresponding horizontal spatial resolution $\Delta s$, as defined in Figure 9, varies with the horizontal range $\varrho$ according to the formula:

$$\Delta s=2\varrho sin\left({\displaystyle \frac{\Delta \phi }{2}}\right)$$

(21)

For an azimuthal resolution of $\Delta \phi =0.{01}^{\circ }$, as used in the experimental part (Section 4), the spatial resolution $\Delta s$ takes the value 0.1745 m at an horizontal distance $\varrho$ of 1 km, 1,745 m at 10 km, and 17.45 m at 100 km. At a distance of about 10 km, the spatial resolution $\Delta s$ is hence about the same as the 2 meters grid resolution of the used DEMs. Figure 10(a) shows the situation at a distance significantly lower than the balance distance of equal resolutions, while Figure 10(b) shows the situation at a distance significantly higher than the balance distance.

In the situation of Figure 10(a), the interpolation ensures that the LAEM gets a value at all discrete azimuth values, also at those situated between DEM points. In the situation of Figure 10(b), several neighbouring DEM points are mapped to the same j value. The condition of listing 1 line 11 however selects only the nearest one. This performs implicitly a subsampling of the DEM points, but without any prior aggregation of those points into a representation of the group. In digital signal processing [55] and image processing [9], such a subsampling is known as performing the unwanted effect of aliasing and has to be preceded by a low-pass filtering to avoid it. Here, we only advice to perform, where it is needed, a correct downsampling of the DEM as preprocessing. However, so far, we have not worked out a solution where a correct downsampling is integrated in the LAEM computation presented in Section 3.4.4.

4. Implementation and Experimental Results

4.1. The swissALTI3D DEM

swissALTI3D [6] is a raster digital Elevation Model (DEM) available currently in the Swiss LV95 [56] (for grid) and LN02 [57] (for height) coordinate systems. All coordinates are in meters. The data exist at two grid resolutions: 0.5 and 2 meters.

The DEM is freely downloadable from [6]. It is split into files (tiles) corresponding to a terrain square of 1 km by 1 km, available in 2 formats: Cloud Optimized GeoTIFF (COG) and zipped ASCII X,Y,Z. ESRI ASCII GRID is also available on request. The selection of the tiles occurs by clicking on a zoomable map of Switzerland (default). Alternative solutions are rectangle and polygonal selection, as well as selection by political regions: canton or municipality. Download of the whole data set is also an option. For a small number of selected tiles (e.g. square of 6 by 6 tiles), the files can be downloaded individually. For larger selections, only a csv file with the list of the file links, e.g. https://data.geo.admin.ch/ch.swisstopo.swissalti3d/swissalti3d_2019_2594-1119/swissalti3d_2019_2594-1119_2_2056_5728.tif, is directly downloadable from the web page.

More information on the DEM is available in a technical report existing in German [58] and in French [59]

4.2. Camera Set-Up and Images

The camera used for the experimentation is an AXIS M5525-E PTZ Dome Network Camera, featuring ${360}^{\circ }$ endless pan, 0-90° tilt, and 10x zoom capabilities. For other camera features, we let the reader refer to the product support web page [60]. The camera has been installed on the roof of the building ENP23 of the HES-SO Valais-Wallis in the town of Sion. This high position gives a ${360}^{\circ }$ view above the neighbouring buildings on the side slopes of the Rhône Valley.

While dome cameras are often installed at a room ceiling to monitor the half-space below the ceiling. In this set-up, the camera has been installed upside down to monitor the half-space above the building roof. To allow a view above surrounding objects on the roof, the camera has been fixed at the top of a 2000 mm high by 800 mm wide frame as shown in Figure 11(a). The frame has been levelled with a spirit level. This mechanical structure has been selected to easily allow, on demand, the installation of other monitoring sensors. As visible in Figure 11(b), a crown-shaped hardware mask has been 3D printed and glued on the camera casing, to avoid collecting privacy sensitive data in taking pictures of the neighbouring buildings.

Figure 12 shows a typical image from the camera. The black area in the bottom of the image is due to the privacy mask. Pictures have been systematically taken with this camera over a whole year at a mean rate of two sessions per day. A subset of the sessions is publicly available from the B2SHARE repository [61] under license “Creative Commons Attribution 4.0 International”. For demanding access to the whole dataset, use the contact information given at [61].

4.3. Determining the Camera Position

As a first step, the camera position has been determined on the online digital map of Switzerland of the Swiss Federal Office swisstopo [62]. Selecting the Aerial imagery as background, it is easy to identify the camera location on the building roof and let the position coordinates being displayed in clicking on the right mouse button. Figure 13 shows the displayed coordinates.

The given elevation value, in the LN02 coordinate system, is of course not the altitude of the camera, but the one of the ground level. To determine the height of the camera position with respect to the ground level, we have used a laser distance meter Leica Disto X4 and found the LN02 altitude of the camera to be 517.5m.

4.4. Preselection and Download of the DEM tiles

We have opted for the COG file format and the 2 meters resolution, since we have estimated this one quite sufficient for modelling the scene monitored by the camera.

Because the part of the Rhône Valley where the camera is located is also a part of the Canton Valais, a simple preselection of the tiles was made by choosing the option “Selection by canton” and then entering the name “Valais”. This generated a csv file to be downloaded. We downloaded the file, which contained, in a single column, a list of 5562 hyperlinks to the individual tile files. The files themselves were downloaded with the help of a MATLAB^® m-file reading the csv file, performing an http request for each hyperlink in the list, and saving the received data as a file into a local directory.

Even if it is not significant, we want to mention that the downloaded DEM files do not fully cover the landscape visible from the camera, since a very small portion of the visible terrain is in France and not surveyed by the Swiss Office of Topography. The concerned region is in a very small angle of view and at more than 40 km from the camera.

4.5. LAEM Computation

4.5.1. Computation of the Geodetic Coordinates

As mentioned in Section 4.1, the swissALTI3D DEM uses Swiss specific coordinate systems for both the horizontal and the vertical datum. In Section 4.3, the camera position has been determined in the same coordinate systems. The transformation into the local AER coordinate system described in Section 3.2 being from a geodetic coordinate system, we have, for using it, first to transform the DEM and the camera position into a geodetic coordinate system.

The LN02 elevation values are orthometric heights. Their conversion into geodetic heights cannot be performed by a mathematical formula, because they are of different natures. The first elevation is linked to gravity, while the second one is linked to geometry. The conversion into the geodetic heights necessitates the official Swiss geoid model CHGeo2004 [63]. The Swiss Federal Office of Topography swisstopo provides the free online service REFRAME to perform coordinate transforms, in particular the transformation of the LN02 heights into the geodetic heights on the Bessel 1841 ellipsoid. This is the transformation we need.

Swisstopo offers different possibilities for using the REFRAME service. First, the transformation can be performed via the web page [64], where a file can be uploaded, the transformation selected, and the resulting file downloaded. This is an easy way for transforming one or a few files. It is also possible to perform coordinate transforms via the REST API [65], but one http request has to be done for each pair or triplet of coordinates. Swisstopo also provides software libraries (DLL, JAR) [66] to enable the use of the REFRAME service within software applications. However, the provided functions also perform the coordinate transform for a single point only. The last offered possibility is the use of the GeoSuite software for Windows [67], provided freely by swisstopo. This software allows to perform transforms of multiple files in batch. The module REFRAME and the type of transformation can be easily selected as represented in Figure 14. The only drawback is that is does not accept swissALTI3D files in GeoTIFF format for the needed altitude transform.

For the transformation of the camera altitude from LN02 to geodetic, we used the REFRAME REST API, performing the http request:

      http://geodesy.geo.admin.ch/reframe/ln02tobessel

?easting=2594281.5&northing=1119555.5&altitude=517.5

We received back the following data:

      {"type": "Point", "coordinates": [2594281.5, 1119555.5, 518.1399090414052]}.

The camera altitude on the Bessel 1841 ellipsoid has been then recorded as 518.14 m.

For the altitude transformation of the swissALTI3D files, we chose the batch ability of the GeoSuite software, combined with a transformation of the swissALTI3D files from GeoTIFF to XYZ, before the altitude transformation, and from XYZ to GeoTIFF after the transformation. We made a test with a batch of 500 files and noticed that the transformation was very resource demanding. Since the difference between the LN02 height and the geodetic height is in the magnitude of not more than 10 meters, we have decided to start the processing with the LN02 heights and perform the transformation only for the tiles that actually remain after the selection described in Section 4.5.2 and Section 4.5.3.

For the transformation of the LV95 coordinates into the geodetic latitude and longitude on the Bessel 1841 ellipsoid, for both the camera position and the DEM files, we simply used the projinv function [68] of the MATLAB^® Mapping Toolbox. The projected Coordinate Reference System (CRS) object to be passed to the function as first parameter can simply be retrieved from any swissALTI3D COG file when reading the file with the function readgeoraster [69]. If R is the variable receiving the second output argument of readgeoraster, the value to be passed to projinv is simply R.ProjectedCRS.

4.5.2. Tile-Level Processing

The preselection of the swissALTI3D tiles performed in Section 4.4 is quite broad. Many of the preselected tiles model landscape which is not at all visible from the camera position. To avoid useless processing, a processing of lower computational complexity is performed at the tile level. This processing extracts information with which it will be possible, while filling the LAEM tile by tile, to decide if a new tile has to be processed in the way given by listing 1 or not at all, because it contains no point visible from the camera position.

In Section 4.5.1, we have explained that the transformation of the LN02 heights into geodetic heights by the swisstopo REFRAME service was too much resource demanding to perform it on the whole preselected files. Hence, the processing at tile level is performed with tiles in LN02 heights. If the transform would have been made, the same processing step would have been performed, but with tiles in geodetic heights. We therefore consider both cases.

To avoid confusion between the local ENU and the LV95 coordinate systems, we will use in the following the symbol u for the easting LV95 coordinate (instead of x) and v for the LV95 northing coordinate (instead of y). A swissALTI3D tile is rectangular in the LV95 coordinate system. Let $u_{\min }$, $u_{\max }$ and $v_{\min }$, $v_{\max }$ be the extremal easting, respectively northing coordinates. For adjacent tiles in east direction, $u_{\min }$ of the eastern tile is strictly larger than $u_{\max }$ of the western tile. The difference is equal to the DEM resolution (2 meters for our choice). The same holds for the northing coordinate v.

The pseudo-code of listing 2 describes the processing performed at the tile level. We add the necessary complementary information, referring to the code by the line number.

Listing 2. Tile-level processing
1	for each tile $T_i$:
2	determine $u_{\min },u_{\max },v_{\min },v_{\max }$
3	if $(u_{\min }\le u_0\le u_{\max })\wedge (v_{\min }\le v_0\le v_{\max })$, current tile is $T_0$
4	for each tile in $\left\{T_i\right\}-\left\{T_0\right\}$:
5	determine $h_{\min },h_{\max }$
6	build B according to equation (23)
7	compute $\phi ,\theta ,d$ for each $(u,v,h)\in B$
8	compute ${\phi }_{\min },{\phi }_{\max },{\theta }_{\max },d_{\min }$ with equation (24)

line 2:

$u_{\min },u_{\max },v_{\min },v_{\max }$ are easy to determine. If the tile format is ASCII XYZ, they can be determined directly from the file name. If the tile format is COG, they can also be found in the file metadata.

line 3: 

If the LV95 coordinates of the camera $(u_0,v_0)$ satisfy this condition, the camera position is within the tile currently considered. This tile is denoted by $T_0$. As explained hereinbefore, adjacent tiles do not have common extremal coordinates, so that the camera position cannot be in more than one tile. In contrary, it can happen that the camera position is between two adjacent tiles. In this case, $\left\{T_0\right\}$ is simply the empty set ∅.

line 5: 

Extremal height values $h_{\min }$ and $h_{\max }$ (either LN02 or geodetic, depending on the file content) can be easily determined by scanning once the file. Performing the processing with original tiles, containing LN02 heights, we increase $h_{\max }$ by a safety margin of e.g. 20 meters, in order to ensure not missing visible points when discarding a tile. Having tiles in geodetic heights, $h_{\max }$ can be kept as it is.

line 6: 

With the extremal coordinate values $u_{\min },u_{\max },v_{\min },v_{\max },h_{\min }$, and $h_{\max }$, we build a tile bounding box, which is a rectangular cuboid in the coordinate systems LV95/height. For the reason that will become clear further, we do not only use the 8 vertices of the rectangular cuboid, but also the edges. The set of $(u,v,h)$ coordinates is hence:

$$\begin{array}{c}B=((\{u_{\min },u_{\min }+\Delta u,\dots ,u_{\max }\}\times \{v_{\min },v_{\max }\})\cup \hfill \\ \hfill (\{u_{\min },u_{\max }\}\times \{v_{\min },v_{\min }+\Delta v,\dots ,v_{\max }\}))\times \{h_{\min },h_{\max }\}.\end{array}$$

(23)

line 7: 

The $(u,v,h)$ coordinates of the points belonging to B are transformed into the values $\phi$, $\theta$, and d of the local coordinate system, where d is the horizontal distance $\varrho$ in the case of the plane approximation of the DEM horizontal datum and d is the slant range r in the case of the spheric approximation. Note that the height is not transformed from LN02 to geodetic. It is kept in the same height coordinate system as the tile.

line 8:line 8:

Following extremal $\phi$, $\theta$, and d values are computed for the points belonging to B:

$$\begin{array}{cc}\hfill {\phi }_{\min }^T& =\underset{(u,v,h)\in B}{min}\phi (u,v,h),\hfill \end{array}$$

(24a)

$$\begin{array}{cc}\hfill {\phi }_{\max }^T& =\underset{(u,v,h)\in B}{max}\phi (u,v,h),\hfill \end{array}$$

(24b)

$$\begin{array}{cc}\hfill {\theta }_{\max }^T& =\underset{(u,v,h)\in B}{max}\theta (u,v,h),\hfill \end{array}$$

(24c)

$$\begin{array}{cc}\hfill d_{\min }^T& =\underset{(u,v,h)\in B}{min}d(u,v,h).\hfill \end{array}$$

(24d)

We use T for tile in the notation in order to not confuse these $\phi$ and $\theta$ values with the extreme LAEM azimuth and elevation values introduced in Section 3.3.3.

Without demonstration, we claim that the values ${\phi }_P$, ${\theta }_P$, and $d_P$ of any tile point P satisfy:

$$({\phi }_{\min }\le {\phi }_P\le {\phi }_{\max })\wedge ({\theta }_P\le {\theta }_{\max })\wedge (d_P\ge d_{\min })$$

(25)

Note that the statement is not true, if B is reduced to the 8 vertices of the rectangular cuboid.

4.5.3. LAEM Computation with Tile Sorting and Filtering

In this section, we describe how the LAEM is computed tile by tile, making use of the information extracted from the tiles, to decide if the current tile has to be processed in the way given by listing 1 or not at all, because it does not contain any point visible from the camera position.

The condition for the decision uses from the tiles the four values computed by Equation (6): ${\phi }_{\min }^T$, ${\phi }_{\max }^T$, ${\theta }_{\max }^T$, and $d_{\min }^T$, and from the LAEM the apparent horizon, i.e. the coordinates of highest elevation for which the LAEM value is not $+\infty$. The azimuth values ${\phi }_{\min }^T$ and ${\phi }_{\max }^T$ determine an azimuthal sector, which we will denote by ${\Phi }^T$. Then, the condition for a tile to be processed is that either the maximal tile elevation ${\theta }_{\max }^T$ is larger than the minimal elevation of the apparent horizon in the sector ${\Phi }^T$ or the minimal tile distance $d_{\min }^T$ is smaller than the maximal LAEM value on the apparent horizon in the sector ${\Phi }^T$. A formal description of the condition will be provided while commenting the pseudo-code of listing 3.

We keep the information on the apparent horizon in a one-dimensional array $\check{\mathbf{I}}$, the length of which is equal to the number $N_{{\phi }^{\prime }}$ of LAEM columns, containing the indexes of the highest discrete elevations. The array element of position j, denoted by ${\check{i}}_j$, is hence formally defined by the following equation:

$${\check{i}}_j=\left\{\begin{array}{cc}min(i\mid d_{i,j}<+\infty ),\hfill & \mathrm{if}\exists i\mid d_{i,j}<+\infty ),\hfill \\ \mathrm{NaN},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(26)

The minimum function is used instead of the maximum, because of the choice made in Section 3.3.3 of indexing the discrete elevations from the highest (top) to the lowest (bottom).

The pseudo-code of listing 3 describes the LAEM computation tile by tile, checking for each tile if they they have to be processed because they may contain visible data. We add the necessary complementary information, referring to the code by the line number.

Listing 3. LAEM computation with tile sorting and filtering
1	$\phantom{\widehat{\check{i}}}$%Initialise $\mathbf{D}$ and $\check{\mathbf{I}}$
2	$\phantom{\widehat{\check{i}}}$$d_{i,j}=+\infty ,\forall i,j$,     ${\check{i}}_j=\mathrm{NaN},\forall j$
3	$\phantom{\widehat{\check{i}}}$update $\mathbf{D}$ and $\check{\mathbf{I}}$ with $T_0$
4	$\phantom{\widehat{\check{i}}}$for each tile in $\left\{T_i\right\}-\left\{T_0\right\}$ sorted ascending according to $d_{\min }^T$:
5	$\phantom{\widehat{\check{i}}}$    compute $j_{\min }^T,j_{\max }^T,i_{\min }^T$
6	compute ${\widehat{\check{i}}}_{j_{\min }^T,j_{\max }^T},{\widehat{d}}_{j_{\min }^T,j_{\max }^T}$
7	if $(i_{\min }^T<{\widehat{\check{i}}}_{j_{\min },j_{\max }})\vee (d_{\min }^T<{\widehat{d}}_{j_{\min },j_{\max }})$:
8	$\phantom{\widehat{\check{i}}}$       update $\mathbf{D}$ and $\check{\mathbf{I}}$ with $T_i$

lines 

3 and 8: The update of $\mathbf{D}$ is performed in the way described in listing 1. The update of $\check{\mathbf{I}}$ is simply performed in inserting the following instruction after the line 11 of listing 1: if $d_{i_P,j_P}==+\infty ,{\check{i}}_{j_P}=i_P$.

line 4: 

Doing the processing without sorting the tiles would give the same final LAEM but would necessitate to process more tiles.

line 5: 

From ${\phi }_{\min }^T$ and ${\phi }_{\max }^T$, the indexes of the corresponding discrete azimuth values $j_{\min }^T$ and $j_{\max }^T$ are computed using the equations (21) and (20a). From ${\theta }_{\max }^T$, the index of the corresponding discrete elevation value $i_{\min }^T$ is computed using the equation (20b).

line 6: 

These are extremal values in the section going from $j_{\min }^T$ to $j_{\max }^T$ of the apparent horizon $\check{\mathbf{I}}$. ${\widehat{\check{i}}}_{j_{\min }^T,j_{\max }^T}$ is the minimal elevation in the section, formally:

$${\widehat{\check{i}}}_{j_{\min }^T,j_{\max }^T}=\left\{\begin{array}{cc}\underset{j_{\min }^T\le j\le j_{\max }^T}{max}({\check{i}}_j\mid {\check{i}}_j\ne \mathrm{NaN}),\hfill & \mathrm{if}\exists j\mid (j_{\min }^T\le j\le j_{\max }^T)\wedge ({\check{i}}_j\ne \mathrm{NaN}),\hfill \\ \mathrm{NaN},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(27)

${\widehat{d}}_{j_{\min }^T,j_{\max }^T}$ is the maximal LAEM value in the section, formally:

$${\widehat{d}}_{j_{\min }^T,j_{\max }^T}=\left\{\begin{array}{cc}\underset{j_{\min }^T\le j\le j_{\max }^T}{max}(d_{{\check{i}}_j,j}\mid {\check{i}}_j\ne \mathrm{NaN}),\hfill & \mathrm{if}\exists j\mid (j_{\min }^T\le j\le j_{\max }^T)\wedge ({\check{i}}_j\ne \mathrm{NaN}),\hfill \\ \mathrm{NaN},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(28)

line 7: 

This is the formal expression of the condition expressed with words at section start.

4.6. Results

The MATLAB^® implementation and the results date from the project time (Dec. 2021-Sept.2023). Only the planar approximation of the DEM horizontal datum has been implemented. The LN02 to geodetic height transformation with REFRAME was tested, also the JAR library, but performed neither on the whole set of the preselected tiles, nor on the set of the finally selected ones.

Figure 15, dating from Sept. 2023, shows the greyscale coded LAEM obtained from the camera position with an angular resolution of $0.{01}^{\circ }$ in azimuth $\Delta \phi$, as well as in elevation $\Delta \theta$. The range in azimuth $[{\phi }_{\min }^{\prime },{\phi }_{\max }^{\prime }[$ is $[0^{\circ },{360}^{\circ }[$ and the one in elevation $[{\theta }_{\min },{\theta }_{\max }]$ is $[0^{\circ },{20}^{\circ }]$. The north-east (azimuth ${60}^{\circ }$) to south-west (azimuth ${240}^{\circ }$) orientation of the Rhône Valley can be easily recognised. The highest elevation of the apparent horizon has a value of about ${18}^{\circ }$ and is southwards, at an azimuth of about ${170}^{\circ }$.

Figure 16 represents in an image the preselected as well as the finally selected tiles. Each tile is represented by a small square. The tile containing the camera position is the isolated white square, close to the image centre. The tiles not containing the camera position, which have been fully processed are in red. There are 805 such tiles. The tiles which have been selected, but finally not processed for filling the LAEM are in black.

5. Discussion

In this paper, I have proposed a new representation of a digital elevation model (DEM) called local azimuth-elevation map (LAEM) and I have described a method to compute this representation. The method has been implemented in MATLAB^® and tested on the swissALTI3D DEM published by the Swiss Federal Office of Topography swisstopo. The implementation has required completing the general method with processing specific to the DEM delivered by swisstopo. That is transforming the Swiss specific coordinates into geodetic coordinates and dealing with the tiling of the DEM. Both, the method principle and its implementation, have been described with sufficient details in order to be analysed on correctness and be exactly reproduced by any interested reader.

The LAEM representation and the method for its computation have been developed in a research project having as goal: the calibration of a PTZ camera with the help of DEM data. Notwithstanding the goal, the LAEM features per se the following advantages:

• It contains only and totally the portion of a DEM which is visible from a given position.
• It is not a subset in any order, but DEM points are arranged with a spatial consistency.
• The dimensionality is reduced from 3D to 2D. The LAEM is comparable to an image of the DEM but with the distance information (horizontal distance $\varrho$ or slant range r) and another 2D coordinate system.
• The scalar value in the 2D space, the slant range r in the spheric approximation case, is rotation invariant.

Regarding the project goal (camera calibration), where the key problem was to find points in the camera image and points in the DEM which correspond to each other, the new representation is an important step on the way to the problem solution for the following reason: it narrows the set of DEM candidates and arranges them in 2D in a way consistent with the arrangement of their projection in the image . For achieving the goal, I launch the following ideas, opening paths for future research:

Depth 

map The LAEM is a depth map, it is similar to the output of a LiDAR. Techniques from the fusion of LiDAR and camera data [70] can be used. In particular, one could work on the similarity between $C^0$ discontinuities [71] of the LAEM and image edges [9]. Note that this is a generalisation of the approach used by Portenier et al. [38], who used the similarity of the apparent horizon in the image and the one in a synthetic image generated from the DEM.

GIS 

augmented LAEM The final goal of the camera calibration is the georeferencing of the camera image to be able to enrich the image with GIS information. Changing the ordering of the procedure, i.e. starting by enriching the LAEM with GIS information, e.g. roads, public lights,…could help in finding correspondences between the camera image and the LAEM. The use of an orthophoto to complement the LAEM, the drawback of which has been discussed in Section 2, would also belong to this approach.

Self- 

calibration In the case of a PTZ camera, a self-calibration can be performed. With it comes along a 3D point cloud. The camera pose is relative to this 3D point cloud. It is then possible to align the 3D point cloud with the LAEM and, by the way, get the pose of the camera in a geographic coordinate system. For aligning the 3D point cloud generated by the self-calibration with the LAEM, the well-known RANSAC algorithm [9] can be used. Another idea is to transform the 3D point cloud in AER coordinates and use the rotation invariance of the range component r and the LAEM value to preselect points which can correspond to each other.

Last but not least, as mentioned in Section 3.4.5, I propose to integrate a downsampling procedure in the method presented for computing the LAEM.