# Terrain descriptors for landscape synthesis, analysis and simulation

O. Argudo, Eric Guérin, H. Schott, Eric Galin

► To cite this version:

O. Argudo, Eric Guérin, H. Schott, Eric Galin. Terrain descriptors for landscape synthesis, analysis and simulation. Computer Graphics Forum, 2025, 44 (2), pp.e70080. 10.1111/cgf.70080. hal-05246144

HAL Id: hal-05246144

https://hal.science/hal-05246144v1

Submitted on 9 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Terrain descriptors for landscape synthesis, analysis and simulation

O. Argudo<fnref n="1" />, E. Guérin<fnref n="2" />, H. Schott<fnref n="3" />, and E. Galin<fnref n="3" />

<figure><img src="image_2.png" /><figcaption>**Figure 1:** A typical example of terrain authored by artists for a AAA game: given an input low-resolution model composed by assembling majors landforms, a high-resolution elevation model augmented with details is generated after computing the slope, sky illumination, accessibility, and drainage, not only to amplify the elevation but also synthesize vegetation occupancy and plausible village settlements.</figcaption></figure>

## Abstract

Synthetic landscape generation is an active research area within Computer Graphics. Algorithms for terrain synthesis and ecosystem simulations often rely on simple descriptors such as slope, light accessibility, and drainage area. Typically, the results are assessed from a perceptual standpoint, focusing primarily on visual plausibility. Other fields, such as Geomorphology and Earth Sciences, have already proposed several analytical descriptors to measure various terrain properties. This work aims to bridge the gap between these disciplines and Computer Graphics. We provide a comprehensive review of commonly used terrain metrics that may be relevant for landscape synthesis, analysis, or simulations. Additionally, we compare the approaches used in Computer Graphics to see if these metrics, or similar ones, have already been introduced. Moreover, we report feedback from a preliminary study conducted with a group of artists to evaluate the potential applications of previously unused metrics. By implementing all these metrics, we enable performance comparisons. Together with the provided correlation matrix, this helps identify instances where a simpler and faster metric can serve as a proxy for a more computationally intensive one.

## CCS Concepts

* **Computing methodologies** → **Shape modeling**;

## 1. Introduction

Terrain synthesis has been an active area of research in Computer Graphics. Although the first approaches relied on procedurally defined fractal noise functions that resembled peaks and valleys [KMN88, MKM89], the emphasis progressively turned towards more authentic modeling either by simulating the underlying geological processes [KMN88, RPP93, NWD05] or by imitating or assembling landforms structures from exemplars gathered from real elevation data [ZSTR07, GDG*17]. A complete description of the vast variety of existing terrain generation methods is beyond the scope of this paper and a state-of-the-art can be found in [GGP*19]. Recent results obtained by Computer Graphics re-

searchers demonstrate that visually appealing, hydrologically coherent, and sometimes geomorphologically convincing scenic landscapes can be synthesized and authored according to artists’ intent.

Nonetheless, evaluating the output of terrain synthesis methods remains a challenging problem. Beautiful rendering of landscapes is sometimes complemented with a perceptual study validating that a proposed method seems more plausible than previous competitors. Nevertheless, casual observers may not be qualified to evaluate the overall realism of terrain morphology, and it is a nearly impossible endeavor to obtain a sufficiently large group of geologists or geomorphologists to get statistically significant results.

The rationale of this work comes from the observation that Com-

<footnote n="1">¹ Universitat Politècnica de Catalunya, Spain</footnote>

<footnote n="2">² LIRIS, INSA-Lyon, CNRS, France</footnote>

<footnote n="3">³ LIRIS, Université Lyon 1, CNRS, France</footnote>
puter Graphics synthesis approaches often rely on relatively simple metrics, such as gradient, slope, illumination accessibility (see Figure 1), or recently on drainage area computation algorithms, even though a variety of metrics exist in geomorphology for analyzing the characteristics of terrains.

This survey aims to bridge the gap between geomorphology literature and Computer Graphics by providing a comprehensive review of commonly used metrics in geological and Earth sciences that could be useful for terrain analysis and synthesis in CG. Our contributions are as follows: 1) We present a state-of-the-art overview of metrics used for analyzing terrains in geomorphology, organized into different categories: fundamental mathematical operators, visibility, landform analysis, hydrology, orometry, and river/ridge network characterization. 2) For every category, we compare the metrics employed in geomorphology with methods used in Computer Graphics and report artists' feedback about their potential utility; 3) We provide the reference implementation used to obtain the figures and timings. Tables 1 and 2 summarize all the metrics along with their properties and performance, and Figure 52 displays their correlation matrix. In the conclusion, we analyze and discuss how some geomorphological metrics could be effectively used in Computer Graphics.

## 2. Local surface properties

Let us represent a terrain surface $\mathcal{T}$ as a height field function $h(\mathbf{p})$: $\mathbb{R}^2 \to \mathbb{R}$ for every point $\mathbf{p} = (x, y)$ in a spatial domain $\Omega \subset \mathbb{R}^2$. We denote the directional derivatives as $h_x = \partial h/\partial x$, $h_{xy} = \partial^2 h/\partial x \partial y$, $h_{xx} = \partial^2 h/\partial x^2$, and so on. The surface area of $\Omega$, expressed in m² or km² depending on the size of the domain, is denoted as $A_\Omega$, whereas the area of the surface of the elevated terrain will be denoted as $A_\mathcal{T}$. We denote $\max_{\mathbf{p} \in \Omega}(\mathbf{h}) = \max_{\mathbf{p} \in \Omega}(\mathbf{p})$, and use similar notations for the minimum. $\tilde{h}(\Omega)$ denotes the mean value of elevations inside $\Omega$, and $\Delta h(\Omega) = \max h(\Omega) - \min h(\Omega)$ the range of elevation within $\Omega$.

Digital Elevation Models (DEM), often referred to as discrete heightfields or elevation maps in Computer Graphics, are the most common representation for terrains. They consist of a set of altitudes arranged on a regular 2D grid. A regular square grid is defined by its rectangular region $\Omega$ with $n$ subdivisions. In this representation, only the altitudes $z_{ij} = h(p_{ij})$ are stored at the vertexes of the grid $p_{ij}$ with the size and resolution of the domain $\Omega$. The distance between point samples will be referred to as $\delta > 0$. While the heightfield representation cannot model overhangs, arches or caves [PGP*19], it is the most widely used nowadays as cartography and topography data is provided in the form of DEM.

Throughout this paper, the properties are expressed whenever possible in a continuous fashion, i.e. we can compute the metric at any position $\mathbf{p} \in \Omega$ and obtain a metric map over the domain $\Omega$. In many cases however, the metrics operate locally over a given neighborhood typically prescribed as a circle of a given radius around the

sample location. Computations are often performed over a discrete squared moving window $\Omega_n$ of size $(2n+1) \times (2n+1)$ centered at $\mathbf{p}$ (Figure 2). Thus $\Omega_1$ represents the classical 3 × 3 window. Here, we omit $\Omega_n$ dependency from the variable names and function parameters whenever evident out of clarity.

### 2.1. Gradient, slope and orientation

The gradient computed from the first order derivatives contains information about the orientation and steepness of the terrain. The gradient of a height field at point $\mathbf{p}$ indicates the steepest direction and is defined as:

$$ \nabla h = \left( \frac{\partial h}{\partial x}, \frac{\partial h}{\partial y} \right) = (h_x, h_y) $$

It is usually approximated using finite differences.

The norm of the gradient corresponds to the slope of the terrain, which can also be expressed as an angle $\beta$:

$$ s(\mathbf{p}) = \|\nabla h(\mathbf{p})\| = \sqrt{h_x^2 + h_y^2} \quad \beta(\mathbf{p}) = \arctan s(\mathbf{p}) $$

Note that the directional slope, i.e., the slope between two points $\mathbf{p}$ and $\mathbf{q}$ is defined as:

$$ s(\mathbf{p}, \mathbf{q}) = \frac{|h(\mathbf{p}) - h(\mathbf{q})|}{\|\mathbf{p} - \mathbf{q}\|} $$

The average slope at a point slightly differs from the slope $s$ and is defined as the average of directional slopes in many directions around $\mathbf{p}$:

$$ \bar{s}(\mathbf{p}) = \frac{1}{n} \sum_{i=1}^{n} s(\mathbf{p}, \mathbf{p} + \epsilon \mathbf{d}_i) \quad \mathbf{d}_i = (\cos 2\pi i/n, \sin 2\pi i/n) $$

In many Computer Graphics applications, we need to compute the surface normal at $\mathbf{p}$:

$$ \mathbf{n}(\mathbf{p}) = \frac{(-\nabla h, 1)}{\sqrt{1 + (\nabla h)^2}} = \frac{(-h_x, -h_y, 1)}{\sqrt{1 + h_x^2 + h_y^2}} $$

The terrain aspect (Figure 3) is the direction of steepest slope opposite to the direction of the gradient. It is often given as the azimuth, which in cartography is defined as the clockwise angle with respect to the North direction. Assuming the $\mathbf{y}$ and $\mathbf{x}$ axes point

<figure><img src="image_3.png" /><figcaption>Figure 2: Window $\Omega_n$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Slope and aspect.</figcaption></figure>
towards the North and the East respectively, the aspect angle is defined as [GW00]:

$$ \alpha(\mathbf{p}) = \pi - \arctan\left(\frac{h_y}{h_x}\right) + \frac{\pi}{2} \frac{h_x}{|h_x|} = \operatorname{atan2}(h_x, h_y) $$

Note that y and x axes are swapped with respect to the mathematical definition of atan2. Since $\alpha(\mathbf{p})$ exhibits a discontinuity between 0 and $2\pi$, some authors [WOB*07] prefer to use the notions of eastness and northness as the two components that define orientation continuously: $\sin(\alpha(\mathbf{p}))$ and $\cos(\alpha(\mathbf{p}))$, respectively.

## 2.2. Laplacian

Second derivatives, or the rate of change of slopes, are widely used in terrain analysis. The Laplacian is defined as the divergence of the gradient, or the sum of unmixed second derivatives:

$$ \Delta h = \nabla \cdot \nabla h = h_{xx} + h_{yy} $$

The Laplacian $\Delta h(\mathbf{p})$ measures how much the average of $h$ in a small sphere centered at $\mathbf{p}$ deviates from $h(\mathbf{p})$, so it can be interpreted as a measure of curvature (see Section 2.3 for a discussion about curvatures).

<figure><img src="image_5.png" /><figcaption>Figure 4: Laplacian and fractional Laplacian.</figcaption></figure>

Unlike the Laplacian, which is a local operator, the fractional Laplacian (Figure 4) represents a key example of the broad category of non-local operators. It can be interpreted as measuring a weighted average rate of change of $h$ across all scales, incorporating contributions from distant points and reflecting long-range dependencies. There are many equivalent definitions of the fractional Laplacian in $\mathbb{R}^n$, we refer the reader to [DL22] for a short survey. Here we use the secondly most used definition in the literature in the form of a well-defined integral on $\mathbb{R}^2$ as follows:

$$ (-\Delta)^s h(\mathbf{p}) = -\frac{1}{2} C(n, s) \int_{\mathbb{R}^n} \frac{h(\mathbf{p}+\mathbf{y}) - 2h(\mathbf{p}) + h(\mathbf{p}-\mathbf{y})}{\|\mathbf{y}\|^{n+2s}} d\mathbf{y} $$

Since the function being integrated is symmetric, it is possible to optimize the computation by evaluating the integral on half of the domain. Moreover, the normalizing constant $C(n, s)$ may be omitted, as many applications use normalized fractional Laplacian maps. Please refer to Appendix A for additional details about $C$ and the relationship with the usual Laplacian operator.

## 2.3. Curvatures

The curvature at a point $\mathbf{p}$ on a curve corresponds to the inverse of the radius of the osculating circle, the circle that better approximates the curve at that point. On a surface, however, the curvature

<figure><img src="image_6.png" /><figcaption>Figure 5: Example of curvatures computed on a height field using fitted interpolating quadrics on $\Omega_3$.</figcaption></figure>

depends on the direction chosen, since there are infinitely many curves that pass through a point. Figure 5 illustrates some measures on a terrain and shows how curvature can be used to identify ridges and the bottom of some valleys.

<figure><img src="image_7.png" /><figcaption>Figure 6: Illustration of different types of curvatures.</figcaption></figure>

Given a normal cross-section at a point $\mathbf{p}$, the curve defined by the intersection between the surface and a plane containing $\mathbf{p}$ and $\mathbf{n}(\mathbf{p})$, we can measure a normal curvature $\kappa_n$. The two principal curvatures correspond to the maximum and minimum values of all the possible normal curvatures: $\kappa_{max}$ and $\kappa_{min}$. Their product defines the Gaussian curvature:

$$ \kappa_G = \kappa_{\max} \cdot \kappa_{\min} = \frac{h_{xx}h_{yy} - h_{xy}^2}{(h_x^2 + h_y^2 + 1)^2} $$

and their arithmetic mean as the mean curvature:

$$ \kappa_M = \frac{1}{2}(\kappa_{\max} + \kappa_{\min}) = \frac{(1+h_y^2)h_{xx} - 2h_{xy}h_xh_y + (1+h_x^2)h_{yy}}{2(h_x^2 + h_y^2 + 1)^{3/2}} $$

From these, we can derive the maximal and minimal curvatures:

$$ \kappa_{\min} = \kappa_M - \sqrt{\kappa_M^2 - \kappa_G} \qquad \kappa_{\max} = \kappa_M + \sqrt{\kappa_M^2 - \kappa_G} $$

Shary [Sha95] defined some curvatures that depend on the slope direction, i.e., that represent features affected by gravitational processes (Figure 6). Profile curvature, also called vertical curvature, is the curvature of the line created by the intersection of a vertical plane in the slope direction. It measures the rate of change of slope,
i.e., the acceleration and deceleration of flow.

$$ \kappa_{\text{pro}} = - \frac{h_{xx}h_x^2 + 2h_{xy}h_xh_y + h_{yy}h_y^2}{(h_x^2 + h_y^2)(h_x^2 + h_y^2 + 1)^{3/2}} $$

Contour curvature, also named plan curvature, measures the curvature of the line created by the intersection of the horizontal plane and the surface. It provides information about the convergence or divergence of flow.

$$ \kappa_{\text{con}} = - \frac{h_{xx}h_y^2 - 2h_{xy}h_xh_y + h_{yy}h_x^2}{(h_x^2 + h_y^2)^{3/2}} $$

Tangential curvature, also called horizontal curvature, is a normal curvature measured from the intersection curve between the surface and the normal plane containing a direction perpendicular to the slope aspect. As reported by [GW00]: it is more appropriate to studying flow convergence/divergence because it does not take extremely large values when slope is small.

$$ \kappa_{\text{tan}} = - \frac{h_{xx}h_y^2 - 2h_{xy}h_xh_y + h_{yy}h_x^2}{(h_x^2 + h_y^2)(h_x^2 + h_y^2 + 1)^{1/2}} $$

There are multiple other curvatures that can be found in the literature, such as rotor, unsphericity, ring curvature, or horizontal/vertical excess curvature. We refer the reader to [SSM02, SEB03, Flo16] for their definition and formulation, and [SEB03] for a comparison of interpolation methods used to compute them from discrete grids.

## 2.4. Quadric surface fitting

Some techniques to compute derivatives and curvatures involve fitting a possibly low degree polynomial surface near p. Evans [Eva80] proposed a variate quadratic surface: $z = a_{20}x^2 + a_{02}y^2 + a_{11}xy + a_{10}x + a_{01}y + a_{00}$

The parameters can be estimated using least squares. For a 3 × 3 neighbourhood and square grid distance between points δ, the exact expression is a function of the nine cell values and their spatial scale (Figure 7).

<figure><img src="image_8.png" /><figcaption>Figure 7: Notations for quadric fitting</figcaption></figure>

$$ a_{20} = ((h_1 + h_3 + h_4 + h_6 + h_7 + h_9) - 2(h_2 + h_5 + h_8))/6\delta^2 $$

$$ a_{02} = ((h_1 + h_2 + h_3 + h_7 + h_8 + h_9) - 2(h_4 + h_5 + h_6))/6\delta^2 $$

$$ a_{11} = (h_3 + h_7 - h_1 - h_9)/4\delta^2 $$

$$ a_{10} = (h_3 + h_6 + h_9 - h_1 - h_4 - h_7)/6\delta $$

$$ a_{01} = (h_1 + h_2 + h_3 - h_7 - h_8 - h_9)/6\delta $$

$$ a_{00} = (2(h_2 + h_4 + h_6 + h_8) - (h_1 + h_3 + h_7 + h_9) + 5h_5)/9 $$

One way to vary the spatial scale of the fitted surface is to re-sample the terrain such that the given 3 × 3 window Ω₁ has the desired scale. Instead, [Woo96] proposed a solution that generalizes to n × n windows and fits a quadric to the original data on a larger extent (see Appendix B for details).

## 2.5. Relief, ruggedness and roughness

The norm of the gradient or slope is not the only metric that is used to model steepness of a terrain (see section 3.2 in [GW00]). A direct metric is the elevation range ΔhΩₙ inside a square domain Ωₙ, sometimes referred to as local relief:

$$ \Delta h_{\Omega_n}(\mathbf{p}) = \max_{\mathbf{q} \in \Omega_n} h(\mathbf{q}) - \min_{\mathbf{q} \in \Omega_n} h(\mathbf{q}) $$

We can define the local mean and the local variance of the elevation as, respectively:

$$ \bar{h}(\mathbf{p}) = \frac{1}{(2n+1)^2} \sum_{\mathbf{q} \in \Omega_n} h(\mathbf{q}) $$

$$ \sigma_z^2(\mathbf{p}) = \frac{1}{(2n+1)^2 - 1} \sum_{\mathbf{q} \in \Omega_n} (h(\mathbf{q}) - \bar{h}(\mathbf{p}))^2 $$

The commonly used Topographic Position Index (TPI) [Wei01] (see Figure 8) is defined as the difference between elevation and the mean elevation [GW00]:

$$ t(\mathbf{p}) = h(\mathbf{p}) - \bar{h}(\mathbf{p}) $$

Since its range of values depends on the window and the terrain, sometimes it is normalized using the standard deviation so that most of its range fits in the [-1, 1] interval:

$$ \hat{t}(\mathbf{p}) = \frac{t(\mathbf{p})}{\sigma_z(\mathbf{p})} = \frac{h(\mathbf{p}) - \bar{h}(\mathbf{p})}{\sigma_z(\mathbf{p})} $$

<figure><img src="image_9.png" /><figcaption>Figure 8: TPI computed with 1500 and 6000 m radius.</figcaption></figure>

Riley et al. [RDE99] introduced the Terrain Ruggedness Index (TRI):

$$ r(\mathbf{p}) = \left( \sum_{\mathbf{q} \in \Omega_n} (h(\mathbf{q}) - h(\mathbf{p}))^2 \right)^{0.5} $$

Although the formula looks very similar to the standard deviation, the Ruggedness Index uses h(p) instead of h̄(p). It is possible to define this index using absolute values instead of squares [WOB*07] over Ωₙ* = Ωₙ - p:

$$ \bar{r}(\mathbf{p}) = \frac{1}{n^2} \sum_{\mathbf{q} \in \Omega_n^*} |h(\mathbf{q}) - h(\mathbf{p})| $$

Jenness et al. [Jen04] defines rugosity (Figure 9) as the ratio between the surface area Aτ of the terrain in the squared window Ω₁
and its projected planar area $A_Ω$:

$$\rho(\mathbf{p}) = \frac{A_T(\mathbf{p})}{A_Ω(\mathbf{p})}$$

Computing the surface area requires building the triangle mesh connecting the centers of each cell with its 8 neighbors and then summing the area of every triangle, whereas the planar area is directly obtained from the extent of $Ω_1$.

<figure><img src="image_10.png" /><figcaption>Figure 9: Rugosity.</figcaption></figure>

Lindsay et al. [LNF19] differentiate between the variability in the elevations — which they refer to as *ruggedness*, in a similar sense as the TRI and local relief metrics described above — and *surface roughness*, which measures the complexity of the surface. These two properties are complementary: we can have a relatively flat but complex terrain (low ruggedness, high roughness) or the opposite, a smooth steep slope (high ruggedness, low roughness). Therefore, they define roughness at $\mathbf{p}$ as the dispersion of surface normal directions in the neighborhood $Ω_n$, and measure it as the spherical standard deviation:

$$\sigma_s(\mathbf{p}) = \sqrt{-2 \ln \|\tilde{\mathbf{n}}\|} \cdot \frac{180}{\pi} \qquad \tilde{\mathbf{n}} = \frac{1}{(2n+1)^2} \sum_{\mathbf{q} \in \Omega_n} \mathbf{n}(\mathbf{q})$$

The vector $\tilde{\mathbf{n}}$ is the average of the $(2n+1)^2$ unit normals in $Ω_n$. The roughness is equal to 0 for flat or inclined planes in which all normals point to the same direction, and increase with surface texture complexity and normals divergence.

As pointed out in [Smi14, LNF19], the term *roughness* is often used as a generic term referring to properties related to local relief or ruggedness, to properties measuring surface complexity, and even combinations of both factors. Figure 10 illustrates the metrics described in this section. For a comparison of additional ruggedness and roughness metrics, we refer the reader to [WYL18].

<figure><img src="image_11.png" /><figcaption>Figure 10: Different relief descriptors computed on $Ω_{16}$.</figcaption></figure>

## 2.6. Asymmetry

By combining the aspect with other metrics, asymmetries can be studied. For example, [PPFB12] use the slope to compute *Hillslope Asymmetry*. Given an aspect direction $θ$, *Hillslope Asymmetry* is defined as:

$$HA(\mathbf{p}) = \log_{10} \left( \frac{\bar{s}(\mathbf{p}, \theta)}{\bar{s}(\mathbf{p}, \theta + \pi)} \right)$$

where $\bar{s}(\mathbf{p}, \theta) = \text{median}(s(\mathbf{q}))$ for locations $\mathbf{q}$ such that $|\alpha(\mathbf{q}) - \theta| < \theta_0$ and $d(\mathbf{p}, \mathbf{q}) < r$. The parameter $\theta_0$ defines the angular spread around $\theta$ and $r$ the maximum distance or window size around $\mathbf{p}$ (Figure 11). From this definition, it follows that the asymmetry of other continuous metrics correlated to orientation can be analyzed.

<figure><img src="image_12.png" /><figcaption>Figure 11: North-South and East-West hillslope asymmetry for $r = 5$ km, scale goes from positive to negative, white is symmetric.</figcaption></figure>

In other cases, asymmetry has been studied with respect to a given linear landform feature, such as a river path or a ridge curve, which requires delineating this features as described in Section 7.1. An asymmetry index for a given property $\mathcal{X}$ can be defined as suggested in [SS20, JRS*22]:

$$Asymmetry(\mathcal{X}) = \left| \frac{\mathcal{X}_L - \mathcal{X}_R}{\mathcal{X}_L + \mathcal{X}_R} \right|$$

where $\mathcal{X}_L$ and $\mathcal{X}_R$ are the values of the property measured at either side (left/right) of the feature.

## 2.7. Analysis

Simple metrics such as gradient and slope are widely used in Computer Graphics in the context of terrain generation, ecosystem simulation, virtual city and road network generation.

The slope is at the heart of most erosion techniques such as hydraulic and thermal erosion, we refer the reader to the state-of-the-art review by Galin et al. [GGP*19] for terrain generation techniques. Hydraulic erosion usually refers to a set of erosive processes involving the transport and deposition of sediments, whereas thermal erosion often refers to a group of various erosional processes other than hydraulic erosion. Both have been widely used since the pioneer works [MKM89, RPP93]. The Stream Power Equation combines the slope and the drainage area (see Section 5) to compute the amount of erosion [CBC*16, CCB*18]. The hill-slope erosion process tends to transport sediments in the direction of the gradient and can finally be expressed as a diffusion-like process involving the Laplacian [CEG*18]. Thermal weathering encompasses different phenomena but can also be described
as a diffusion process that smoothes landforms, which is often achieved by introducing a Laplacian-based term in the erosion equations [CBC*16, SPF*23]. Aeolian erosive processes often rely on the computation of the wind flow whose direction and velocity are approximated by procedural rules from the elevation above the ground and the local slope. This has been used in the context of snow transport modeling [CEG*18] and desert dune simulation [PPG*20].

The slope is also a natural proxy for approximating the amount of soil deposited on the underlying terrain bedrock as it has the advantage of being local and fast to compute. The presence or density of vegetation can be approximated according to slope since vegetation hardly grows on steep slopes. It is thus widely used in procedural modeling together with altitude ranges in order to populate terrains with specific species. On the ecosystems simulation side, abiotic conditions usually encompass simple characteristics of the terrain such as the altitude or the slope [MHS*19, KGG*20].

Modeling in the domain of the gradient has been successfully used to build an interactive terrain editor based on a constrained Poisson reconstruction [GPM*22]. Finally, the slope is also used in the context of village generation [EBP*12] and particularly road path computation [GPMG10], as steep slopes prevent road and building construction.

Slope, curvature [KS14], Laplacian and particularly fractional Laplacian [Bro14] have been used in Cartography for effective relief shading and enhancing topographic details. A complete overview on cartographic shading techniques is beyond the scope of this paper; we refer the reader to a recent article [JHS*21] and to the survey paper [FSH20] for a comparison of shading techniques. In Computer Graphics, exaggerated shading [RBD06] has been investigated for a wider range of objects.

Local surface properties such as slope or Laplacian are commonly used by artists to guide procedural populating of materials and vegetation. Their comments in our artist feedback study (see Section 8) show that curvatures are not yet exploited and could be used to add high-resolution details onto terrains by providing interesting patterns. Another interesting feedback is the potential usage of the fractional Laplacian to a proxy for sediment deposition maps.

### 3. Visibility

Visibility is a particularly interesting descriptor in many applications such as forest management, landscape planning or archaeology [Llo03] and has been extensively studied in Geographic Information Systems (GIS).

#### 3.1. Viewsheds

The viewshed of a point $V(p)$ is defined as the set of points visible from $p$. This definition relies on the line-of-sight between two points: $p$ and $q$ are mutually visible if the segment $pq$ does not have any intersections.

In the case of terrain analysis, the potential occluder is typically the terrain surface. The observer $p$ and target $q$ are often points

on the surface of the terrain $T$. Some applications may consider a point above it, e.g., to account for the height of a standing observer or a watchtower. In this case, visibility may not be reciprocal.

If instead of a single observer we use a set of observer locations $P$, we can define the cumulative viewshed [Whe95] of $P$ as the logical union of individual viewsheds (Figure 12), i.e. the set of points that are visible from at least one location in $P$:

$$V(P) = \bigcup_{\mathbf{p} \in P} V(\mathbf{p}) = \{\mathbf{q} \in T \mid \exists \mathbf{p} \in P, \mathbf{q} \in V(\mathbf{p})\}$$

The visibility index $v_{idx}(p)$ allows to compare different observer locations by measuring the area (or locations count) seen from $p$:

$$v_{idx}(\mathbf{p}) = \sum_{\mathbf{q} \in V(\mathbf{p})} A(\mathbf{q})$$

$A(q)$ represents the area of $q$, or 1 in the case of counting discrete locations. In certain contexts, the reciprocal definitions of viewshed and visibility index are applied, where points $q$ indicate the locations from which $p$ is visible.

Fisher [Fis92] introduced the concept of fuzzy viewshed to account for potential uncertainty in $h$ that may affect visibility. Unlike a binary viewshed, where $q \in V(p)$ strictly denotes that $q$ is visible from $p$, a fuzzy viewshed assigns a value $\tilde{V}(p, q) \in [0, 1]$ to indicate the likelihood of $q$ being seen from $p$. In [Fis92], this likelihood is computed as the proportion of binary viewsheds in which $q \in V(p)$, each viewshed generated from a different perturbed terrain surface $T'$. Viewshed uncertainty can also be applied to the observer, the targets, or atmospheric conditions, for example.

The complexity of a viewshed is intuitively related to the number of local maxima in the domain $\Omega$. However, as demonstrated in [AJL*21], merely counting the peaks is not sufficient, as a small adjustment to the vertical axis followed by a rotation can eliminate all local maxima while still retaining the overall complexity of the viewshed, i.e. viewshed complexity is invariant under affine transformations. Thus, the authors define prickliness as the maximum number of peaks in $\Omega$ across all possible affine transformations.

Since generating viewshed maps can be computationally demanding, advanced algorithms leverage previously computed visibility data from other locations and utilize spatial coherence [TCZR15] to speed up the calculations. Instead of measuring the visible surface, Cervilla et al. [CTV*17] propose an efficient algorithm to compute the volumetric viewshed, i.e., the volume visible from a location.

<figure><img src="image_13.png" /><figcaption>Figure 12: Visibility index of a terrain and viewshed from a point.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: Zenith and Nadir angles $\phi_d(p)$ and $\psi_d(p)$ at a given position $p$ and towards a direction $d$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14: Positive (left) and negative (right) openness.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15: Comparison between accessibility and sky view.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16: Accessibility and sky-view factor computed with 2048 sampling rays per location.</figcaption></figure>

## 3.2. Openness

At a given location, zenith and nadir are defined as the opposite vertical directions that point above and below the terrain respectively in a gravitational sense: i.e. the nadir is aligned with the gravity vector. Consequently, zenith and nadir angles can be defined as angles to these directions.

Yokoyama et al. [YSP02] used these angles to introduce the concept of openness. Given a position $p$ and a horizontal (planar) direction vector $d$, they define the maximum zenith angle $\phi_d$ that a ray can have such that it is tangent to the terrain surface but does not intersect it (Figure 13); the nadir angle $\psi_d$ is defined similarly:

$$ \phi_d(\mathbf{p}) = \frac{\pi}{2} - \max_{\lambda>0} \left\{ \arctan \left( \frac{h(\mathbf{p} + \lambda\mathbf{d}) - h(\mathbf{p})}{\lambda} \right) \right\} $$

$$ \psi_d(\mathbf{p}) = \frac{\pi}{2} + \min_{\lambda>0} \left\{ \arctan \left( \frac{h(\mathbf{p} + \lambda\mathbf{d}) - h(\mathbf{p})}{\lambda} \right) \right\} $$

Positive openness $\Phi(p)$ is the mean of zenith angles in all directions, whereas negative openness $\Psi(p)$ is the mean of nadir angles in all directions (Figure 14). A local implementation limits the ray distance to a user-prescribed maximum for better performance.

## 3.3. Accessibility and Sky-view factor

Openness can be interpreted as the simplification of the obscurance or, inversely, the accessibility: the fraction of a sphere of finite or infinite radius visible from a given location. Alternatively, the Sky-view factor [DF90, ZOK11] considers the visible fraction of sky, i.e. the fraction of unobstructed celestial hemisphere $\mathcal{H}$ above a location (Figure 15). Note that in inclined planes, where an hemisphere oriented with the surface normal would be fully visible, part of the sky is occluded. This definition differs from the typical Ambient Occlusion in Computer Graphics (see Section 3.5 for a discussion).

It is possible to compute the sky-view factor $v_H(p)$ casting rays from $p$ towards directions sampled on the sky hemisphere (Figure 16). Another possibility is integrating the horizon angles around $p$. The horizon angle at $p$ in direction $d$ is defined as the maximum angle toward any other surface in that direction:

$$ \phi_d(\mathbf{p}) = \max_{\mathbf{q}=\mathbf{p}+\mathbf{r}\mathbf{d}, \mathbf{r}>0} \arctan s(\mathbf{p}, \mathbf{q}) $$

Then, integrating over planar directions around $p$ [DF90]:

$$ v_H(\mathbf{p}) = \frac{1}{2\pi} \int_0^{2\pi} \left[ \cos\beta \cos^2\varphi + \sin\beta \cos(\Phi - \alpha) \cdot \left(\frac{\pi}{2} - \varphi - \sin\varphi \cos\varphi\right) \right] d\Phi $$

The angles $\beta$ and $\alpha$ denote the slope and aspect angles at $p$ (see Section 2.1). A simple approximation consists in assuming that the terrain is flat around $p$ with slope angle $\beta(p)$:

$$ \tilde{v}_H(\mathbf{p}) \approx \frac{1 + \cos\beta(\mathbf{p})}{2} $$

## 3.4. Solar irradiance

Surface irradiance can be expressed as [DR95]: $E(p) = E_S(p) + E_A(p) + E_T(p)$, the sum of the direct sun irradiance $E_S$, the diffuse sky hemisphere irradiance $E_A$ and the indirect irradiance from nearby terrain $E_T$ (see Figure 17). Determining actual irradiance involves defining physical properties like solar radiance, atmosphere composition and phenomena such as transmission scattering, which is particularly computationally intensive. An important part of the computation is related to geometric visibility problems driven by terrain topography. For an extended presentation, we refer the reader to [BA09].

Under clear sky, direct irradiance $E_S$ is the dominant term. Sun rays are defined by sun elevation or inclination angle $\theta$ with respect to horizon, and sun azimuth angle $\phi$ with respect to North. The angle between rays with direction $(\theta, \phi)$ and the oriented plane at $p$
<figure><img src="image_18.png" /><figcaption>Figure 17: Terms for solar radiation.</figcaption></figure>

with slope angle $\beta(\mathbf{p})$ and aspect $\alpha(\mathbf{p})$ is:

$$ \cos\gamma = \cos\beta \sin\theta + \sin\beta \cos\theta \cos(\phi - \alpha) $$

Occlusions from nearby relief (shadows) are modeled using the visibility term $\zeta(\mathbf{p}, \omega) \in \{0, 1\}$. Putting it all together, the direct irradiance term is defined as:

$$ E_S(\mathbf{p}) \propto \max(0, \cos\gamma) \cdot \zeta(\mathbf{p}) $$

Diffuse irradiance $E_A(\mathbf{p})$ is caused by atmospheric scattering of incoming rays. Assuming isotropic sky radiance, the sky-view factor serves as good estimator for $E_A(\mathbf{p})$ since it computes the fraction of the hemisphere not occluded by nearby relief:

$$ E_A \propto v_H $$

However, the sky is generally anisotropic - brighter near the horizon and the sun - and accurate solutions might require ray casting the sky hemisphere from $\mathbf{p}$.

Indirect irradiance $E_T$ is the most computationally demanding term. Accurate computation involves evaluating the Bidirectional Reflectance Function of the surface terrain and using Global Illumination algorithms to integrate indirect light reflected from visible terrain. Its contribution to $E$ is typically negligible, except in specific cases such as snow-covered terrain. Under the critical but rarely met assumption of isotropic surrounding relief, the equation simplifies to:

$$ E_T = (\overline{E_S} + \overline{E_A}) \cdot r_0 \cdot v_T $$

where $\overline{E_S}$ and $\overline{E_A}$ are the averaged irradiances over visible terrain, $r_0$ the average terrain reflectivity, and $v_T$ a terrain-view factor estimated as [DF90]:

$$ v_T \approx \frac{1 + \cos\beta}{2} - v_H $$

Finally, even if solar radiation is assumed to be symmetric on eastern and western slopes, the latter are typically warmer due to higher moisture near the ground during morning hours, which increases latent heat flux on east-facing slopes. The Diurnal Anisotropic Heat Index [BA09] approximates this as:

$$ \theta(\mathbf{p}) = \cos(\alpha_{\max} - \alpha(\mathbf{p})) \cdot \arctan\beta(\mathbf{p}) $$

Where $\alpha_{\max}$ represents the direction with maximal heat surplus (e.g., 202.5° for SSW), and $\alpha$ and $\beta$ are aspect and slope angles, respectively (see Figure 18, right). Empirical heat load indices can also be derived from fitted data; for examples, we refer the reader to [MK02].

<figure><img src="image_19.png" /><figcaption>Figure 18: Direct irradiance $E_S$ for latitude 45°N averaged over a year (left), and Diurnal Anisotropic Heat Index (right).</figcaption></figure>

## 3.5. Analysis

Computing the obscurity [SAG98] has been widely studied in Computer Graphics for Global Illumination and Ambient Occlusion [Lan02, MFS09]. Visibility, also referred to as accessibility [Mil94], is frequently used for determining the amount of light accessing a point and shading. Formally, for a point $\mathbf{p}$ with a normal vector $\mathbf{n}$, Ambient Occlusion is defined as the integral over the hemisphere $\mathcal{H}$ of the visibility function $\zeta(\mathbf{p}, \hat{\omega})$ with $\zeta = 1$ if there are no occluders between the point $\mathbf{p}$ and the environment in direction $\hat{\omega}$, and $\zeta = 0$ otherwise:

$$ a(\mathbf{p}, \mathbf{n}) = \frac{1}{\pi} \int_{\mathcal{H}} \zeta(\mathbf{p}, \hat{\omega}) (\mathbf{n} \cdot \hat{\omega}) d\hat{\omega} $$

Accelerated techniques for fast global illumination of terrains received specific attention in [SN08].

In Cartography, [FSH20] presents a survey of relief shading techniques relying on local and global illumination models similarly to Computer Graphics. In particular, clear sky shading [KS14] is a variation of ambient occlusion shading that was preferred over other relief shading methods in a recent user study [FSH20].

Exaggerated shading [RBD06] combined with local light alignment [MVNT21] combine local shading models to enhance features. Similar approaches exist in cartography: a vast amount of literature compares different shading strategies. Instead of using accessibility, they often rely on Laplacian performed at different scales, and on multiple light sources to locally enhance the shading parts of the relief aligned with the direction of the light [KPJ*21].

Direct and indirect illumination impact plant development and are frequently used in ecosystem simulations [CGG*17, MHS*19, GLCC17]. Accessibility and sky-view are also often used as a proxy for evaluating indirect illumination without integrating the trajectory of the sun combined with the average cloud cover over a year [KGG*20]. They are also used for evaluating the influence of lighting over temperature variations depending on the orientation to simulate snow melting [CEG*18]. Large-scale terrain snow cover generation algorithms often proceed by displacing the surface of the terrain or using a layered data-structure to record the thickness of the snow layer, which is derived from the slope, and using local accessibility and occlusion characteristics and snow melting to direct and indirect illumination [PTS99, FB07, RLD15].

Visibility and accessibility is central in analysing and designing a landscape, especially in the video game industry where the user
playability is often linked to what they can actually see. Regarding artist feedback (see Section 8), they mentioned an interesting avenue for openness $\Phi(\mathbf{p})$ to be used as a proxy to snow presence. In addition, they also have detected a potential of viewsheds $V(\mathbf{p})$ to delimit regions, to automatically place villages or castles.

## 4. Landforms

There is no unique definition of a set of basic landform units: they are usually chosen according to the application and detail analysis. A recent review is provided by Mokarram and Sathyamoorthy [MS18]. Landforms are derived from one or several previously presented metrics using segmentation or threshold values.

### 4.1. Curvatures-based

One of the most commonly used landforms set relies on peak, pit, ridge, valley, saddle and plain classes for classifying points $\mathbf{p}$ such that slope is null: $s(\mathbf{p}) = 0$. The classification can be derived from the relative elevations of a $\Omega_1 = 3 \times 3$ neighborhood with respect to the point [PD75]. Wood [Woo96] improved the robustness of this

<figure><img src="image_20.png" /><figcaption>**Figure 19:** Slope units and corresponding landforms.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 20:** Curvature-based landforms from fitted quadrics on $\Omega_{30}$.</figcaption></figure>

classification using curvatures derived from fitted quadrics, larger $\Omega_n$, and tolerance thresholds.

Different combinations of curvatures have been widely used (Figure 20). For example, Dikau *et al.* [Dik89] used plan and profile curvatures to obtain a $3 \times 3$ characterization of slope patterns (Figure 19). For planar areas, the maximum and minimum curvatures correspond to the six classes stated above.

A more complete classification of slope patterns (non-special points) into twelve classes based on the signs of tangential, profile, mean, difference and total Gaussian curvatures is proposed by Shary *et al.* [Sha95, SSM05].

### 4.2. Based on Topographic Position Index

Weiss [Wei01] proposed a classification of landforms based on the *Topographic Position Index* (TPI) at two spatial scales (Figures 21 and 22). To compare the small and large scales, the indices are nor-

<figure><img src="image_22.png" /><figcaption>**Figure 21:** Characterization of landforms using TPI at two scales: the dark lines represent one standard deviation for each scale.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 22:** Landforms from TPI at scales 300 and 2000m.</figcaption></figure>
malized to mean 0 and standard deviation 1. Distinction between flat areas and open slopes (classes 5 and 6 in Figure 21) comes from a threshold on the local slope s, for example 10 degrees.

### 4.3. Ridges and valleys

The percentage of values in a neighborhood Ω that are lower than the central point p yields a simple metric that differentiates pits and valley bottoms with values close to 0 from peaks and ridge tops, with values close to 1 (Figure 23). This metric was called percentile by Gallant and Wilson [GW00], although some studies use the term prominence [Llo01]. To prevent confusion with the prominence of a mountain peak (see section 6), we adopt the term percentile:

$$p(\mathbf{p}) = \frac{|\{\mathbf{q} \in \Omega | h(\mathbf{q}) < h(\mathbf{p})\}|}{|\{\mathbf{q} \in \Omega\}|}$$

<figure><img src="image_24.png" /><figcaption>Figure 23: Percentile metric with medium and large scale radii.</figcaption></figure>

On a larger structural scale, peaks and valleys can be segmented using the White Top Hat and Black Top Hat transformations as proposed by Rodriguez et al. [RMCRD02]. Given a structuring element ε, erosion ⊖ and dilation ⊕ morphological operators are defined as:

$$(h \ominus \varepsilon)(\mathbf{p}) = \inf_{\mathbf{p} \in \Omega_{\varepsilon}} h(\mathbf{p} + \mathbf{q}) \quad (h \oplus \varepsilon)(\mathbf{p}) = \sup_{\mathbf{p} \in \Omega_{\varepsilon}} h(\mathbf{p} + \mathbf{q})$$

Note that in a discrete height field and supporting element ε, inf and sup are the min and max of h in the neighborhood defined by ε: therefore erosion and dilation can be interpreted as a minimum and maximum filter respectively. Recall that the morphological opening and closing operators are defined as: $O^ε = ((h \ominus \varepsilon) \oplus \varepsilon)$ and $C^ε = ((h \oplus \varepsilon) \ominus \varepsilon)$. The White Top Hat W and Black Top Hat B are defined using a user-defined threshold T (see Figures 24 and 25):

$$W^\varepsilon = \{\mathbf{p} | h(\mathbf{p}) - O^\varepsilon(\mathbf{p}) \geq T\} \quad B^\varepsilon = \{\mathbf{p} | C^\varepsilon(\mathbf{p}) - h(\mathbf{p}) \geq T\}$$

### 4.4. Fuzzy and multi-scale classification

The shape and class of a landforms often vary with scale, e.g., a small depression near the top of a mountain peak could be classified either as a pit or peak depending on the analysis scale. For this matter, Fisher and Wood [FW98] proposed using fuzzy memberships instead of reporting a single landform category. The classification is computed following the work of Wood [Woo96] at every scale independently. Then the fuzzy membership to every landform type is defined as the relative number of times it was classified as such. This concept was further developed by Fisher et al. [FWC04]

<figure><img src="image_25.png" /><figcaption>Figure 24: White Top Hat and Black Top Hat transformations.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 25: Black Top Hat ridges (red) and valleys (blue) obtained using Ω7 and two different thresholds.</figcaption></figure>

by defining the fuzzy membership as a weighted average of the individual scale memberships with weights that depend on the scale (although the authors only tried constant weighting). The highest membership defines the dominant or modal class at a given location, and the entropy of the memberships allows for evaluating the corresponding variability.

Schmidt et al. [SH04] apply the concept of fuzzy membership to the curvature thresholds to define the 15 types of form elements as in Figure 19. Larger-scale context is given via the Black Top Hat and White Top Hat transformations [RMCRD02] to classify land as hill, valley, hillslope, or some intermediate fuzzy value if both transformations return non-zero values. Finally, a set of rules combines the landforms with context information to build the final membership maps for land elements (Figure 26): ridge, shoulder, backslope, hollow, spur, terrace, footslope, and valley bottom.

Instead of combining multiple scales, [SA05] proposed multi-scale landform signatures by identifying the relevant ranges in a scaling graph, e.g., the curvature of a fitted polynomial depending on the scale. Ranges are identified from the breakpoints in the derivatives of the scaling graph.

### 4.5. Geomorphons

Introduced by [JS13], Geomorphons are a landform classification system independent of analysis scale (Figure 28). Instead of comparing heights with respect to the central cell to derive different form patterns, they compare the nadir ψd and zenith φd angles along the 8 principal directions [YSP02] (see section 3) to obtain an 8-valued ternary pattern.
<figure><img src="image_27.png" /><figcaption>Figure 26: Landforms following [SH04].</figcaption></figure>

For a given direction d, the pattern Δd in the direction d is computed according to an elevation parameter T as:

$$ \Delta_d = \begin{cases} 1 & \text{if } \psi_d - \phi_d > T \\ 0 & \text{if } |\psi_d - \phi_d| \le T \\ -1 & \text{if } \psi_d - \phi_d < -T \end{cases} $$

For 8 directions, a total of $3^8 = 6561$ ternary patterns are possible (see Figure 27). Taking rotations and symmetries into account reduces this count to 498, and even down to 10 classes (flat, peak, ridge, shoulder, hollow, slope, spur, footslope, valley, and pit) by simply counting the number of 1 and -1 entries in the pattern. Patterns can be efficiently assigned to these meaningful landforms using a look-up table.

<figure><img src="image_28.png" /><figcaption>Figure 27: Example of geomorphon configurations.</figcaption></figure>

The benefit of this method with respect to differential geometry approaches like [Woo96] or [SH04] is that there is no need to apply the analysis at multiple scales. While there is a maximum distance parameter L when computing the zenith and nadir angles, the authors showed results converge to a stable solution with increasing L whereas increasing window size in the other approaches usually changes the identified landforms.

In a follow-up work [JNS14], histograms of geomorphon pairs in neighboring cells were used to define a similarity metric between two terrains. The authors showed applications of it for landscape search, i.e. deriving a similarity map between a large terrain and a query height field, and for segmentation of terrains given a few reference templates per class.

<figure><img src="image_29.png" /><figcaption>Figure 28: Geomorphons computed using maximum distance 1km.</figcaption></figure>

## 4.6. Analysis

Contrary to local surface properties and visibility, landform analysis has been seldom utilized in Computer Graphics. One exception is the use of geomorphons (see Section 4.5) to associate the perceived realism of terrains to the distribution of landforms [RKv*22]. In this work, the authors further employ this information to enhance the authenticity of synthetic terrains.

Note however that the concept of features and patterns has been employed extensively in Computer Graphics without explicitly correlating to the extraction of landforms using traditional geomorphological descriptors. Procedural, erosion simulation and image-based terrain generation models often produce landforms such as ridges [ZSTR07] and riverbeds [PDG*19]. The procedural terrain generation model by Genevaux et al. [GGG*13] relies on the definition of landforms placed and combined in a constructive tree [GGP*15]. Sparse modeling applied to terrains [GDGP16, AAC*17] generate complex elevation by combining relief patches stored in atoms and could benefit from a more semantic representation of the dictionary.

While landforms analysis did not receive much attention in Computer Graphics, it is worth mentioning that the modeling and generation of specific landforms are common applications. Particularly, specific techniques were proposed for generating canyons [MDCPBS14], table mountains [BA05], dunes [PPG*20], and waterfalls [EPCV15].

We collected artist feedback regarding a terrain segmentation based on TPI (Section 4.2). In general, they consider it potentially useful as a selection tool.

## 5. Hydrology

Hydrology studies the movement or distribution of water on a terrain, its processes, and the properties of the flow. A fundamental first step in hydrology analysis is the definition of the flow routing model that describes how water flows to or from a given position
<figure><img src="image_30.png" /><figcaption>Figure 29: Comparison of different flow models from a single source point on a tilted plane.</figcaption></figure>

p (Section 5.1). Hydrology metrics and analysis depend on flow routing, therefore we first review some commonly used flow models. Figure 29 illustrates a synthetic comparison of different flow routing models. Although these apply to regular grid height fields, their concepts can be adapted to other terrain representations, e.g., Triangulated Irregular Networks. We then review depression and flat area detection algorithms (Section 5.2), hydrological metrics (Section 5.3 and 5.4) and erosional indices related to drainage basin dynamics (Section 5.5).

## 5.1. Flow routing models

The deterministic 8-neighbors model [OM84], or simply D8 model, considers the eight neighbors in $\Omega_1$ of a cell p and assigns the flow direction towards the cell with the steepest descent (Figure 30). If all cells are equal or higher than p, no direction is assigned.

One limitation is the occurrence of artificial parallel straight flow lines on regular slopes due to the 45° discretization of directions. The Rho8 model, or stochastic 8-neighbors model [FL91], chooses the outflow cell randomly among the set of lower cells with probabilities proportional to the slope.

While Rho8 solves the straightness of D8, it is not free of objection due to its stochastic nature. The Kinematic Routing Algorithm [Lea92] measures flow from a source as a dimensionless rolling ball following the steepest direction and does not suffer from the quantization artefacts of D8. To compute the flow direction, the elevations at the four corners of a grid cell - averaged from neighbors - are used to compute the best-fit local plane. The outlet is defined as the intersection of the flow direction with the edges of the cell and determines which of the four lower cells the flow propagates to. The outlet position then becomes the flow origin in the downstream cell.

Instead of one fitted plane, the $D\infty$ model [Tar97] determines eight slope directions from the triangles between the cell center and pairs or consecutive neighbors (Figure 29). The flow is then routed towards the steepest direction, but distributed to cells in the corresponding triangle by taking into account the angle between this direction and the cells' centers.

<figure><img src="image_31.png" /><figcaption>Figure 30: Single-direction (D8) and multiple direction (MFD) computation of drainage area A.</figcaption></figure>

Multiple Flow Directions algorithms distribute flow to more than one neighbor. This strategy, also called divergent flow, was proposed to overcome the artefacts of D8 and notably allow diverging flow paths. The model proposed by Holmgren et al. [Hol94] distributes the incoming flow on a cell to all its lower cells $\mathcal{L}(p)$ with proportions $f_i$:

$$f_i = \frac{s(\mathbf{p}, \mathbf{q}_i)^\gamma}{\sum_{j \in \mathcal{L}(\mathbf{p})} s(\mathbf{p}, \mathbf{q}_j)^\gamma}$$

$s(\mathbf{p}, \mathbf{q}_i)$ denotes the slope between p and the i-th neighbor $\mathbf{q}_i \in \mathcal{L}(p)$. This model is a generalization of the one previously proposed by [QBCP91], which corresponds to the case of $\gamma = 1$ and suffers from excessive dispersion of the flow. The larger the value of $\gamma$, the smaller the dispersion of the flow. When $\gamma \to \infty$, the model converges to D8. Holmgren [Hol94] concludes that setting $\gamma \in [4, 6]$ best approximates real flows, while diminishing the risk of outliers.

## 5.2. Depressions and flat areas

A depression is an area of the terrain with no outflow: its closed perimeter is higher than any interior point so the water flow out - if we ignore possible underground drainage systems. While depressions can result from natural features such as pits, dolines, glacial lakes, and endorheic basins, most of the time small depressions appear
<figure><img src="image_32.png" /><figcaption>**Figure 31:** Depression breaching and depression filling on a noise turbulence terrain.</figcaption></figure>

pear in Digital Elevation Models due to errors in acquisition or discretization. The flow algorithms described in Section 5.1 assume strictly decreasing discharges to generate a continuous flow towards the borders of the map or the sea. Depressions or flat areas break this assumption, therefore it is standard to preprocess the height field to remove such unwanted features either by filling or breaching depressions (Figure 31).

Depression filling methods raise the elevation of points in the depression to the elevation of its outlet, i.e., the lowest point in the depression contour. Variations of the Priority-Flood algorithm [WL06, BLM14, ZSF16] employ a priority queue to process cells in increasing elevation order, starting from the borders of the height field. Every time a cell is visited, its unprocessed neighbors are added to the queue. If their elevation is lower than the elevation of the current cell, they are raised to match the same height. To avoid filling depressions with flat areas, a tiny elevation offset $\epsilon$ is added when visiting neighboring depression cells to guarantee a consistent flow routing towards the outlet.

In contrast, depression breaching methods carve a channel from the depression cells towards the outlet or the nearest lower cell. They typically produce smaller changes to the original elevation map, unless the depression is exceptionally deep. The Priority-Flood algorithm can also be adapted to breach depressions [SVC03, BLM14]. When processing neighbor cells, the backwards direction is stored so that once a depression is found, a path can be traced back to an outlet. Then, the elevations of the cells along this path are lowered to the same elevation as the outlet.

Hybrid approaches combine filling and breaching to minimize the alterations on the original height field. The algorithm presented in [Lin16] starts by shallowing single-cell pits, which are likely to be sampling artifacts and cause deep breach incisions, and sets their elevation just slightly below their lowest neighbor. Then, they use a breaching algorithm similar to [SVC03] constrained by a maximum carving depth or breaching channel length. The remaining depressions are eventually filled with a slight slope.

It is also possible to solve depression routing without altering the terrain topography. Cordonnier et al. [CBB19] construct a graph of the basins, with edges connecting adjacent basins through their lowest pass. The Minimum Spanning Tree of this planar graph, which can be computed in linear time, is then used to determine the flow between basins. Recently, this algorithm has been adapted for

GPU computation [JKG*24]. This approach, however, only supports single direction flows.

For the particular case of flat areas with at least one lower outlet, [GM97] resolve drainage direction by combining two gradients: one towards lower terrain, and one away from higher terrain. Similarly, [SVC03] use the inverse distance away from higher cells as a cost weight, and find the minimum cost path from each flat cell towards an outlet.

### 5.3. Drainage basin metrics

Once we have a flow direction for any point $\mathbf{p}$, we can derive several flow metrics that can serve as descriptors or predictors of fluvial and erosional processes.

The most important concept is probably the Drainage Basin $\mathcal{D}(\mathbf{p})$ of a location $\mathbf{p}$, defined as the area of the terrain from where the surface flow can potentially reach $\mathbf{p}$ (Figure 32). The upslope Drainage Area $A(\mathbf{p})$, also called Total Catchment Area or Contributing Area, is defined as the area of $\mathcal{D}(\mathbf{p})$. The drainage basin and total catchment area can be computed by tracing the flow from higher to lower positions following a given propagation model as described in Section 5.1. The choice of flow model can significantly impact the output of drainage area computations, and therefore other metrics derived from it.

<figure><img src="image_33.png" /><figcaption>**Figure 32:** Drainage basin and drainage area.</figcaption></figure>

The Specific Catchment Area $\tilde{A}(\mathbf{p})$ is defined as the upstream catchment area per unit length of contour or cell size. Note that using unit square cells results in $\tilde{A}(\mathbf{p})$ being equal to $A(\mathbf{p})$, but their units and physical interpretation are different.

Flow propagation algorithms can also be used to accumulate upstream metrics in a catchment, for example the sum of upstream elevations or slopes. Dividing these accumulated values by the number of cells yields a basin average, and by further dividing by area cell, a density.

Given delimited drainage basin $\mathcal{D}$ with area $A$, perimeter $P$ and elevation range $\Delta h(\mathcal{D}) = \max h(\mathcal{D}) - \min h(\mathcal{D})$, we can compute different basin shape and relief descriptors. See, for example, Table 1 in [Son17] for a summary and additional references.

Some descriptors reflect the basin shape: Basin Circularity is the ratio between the area of the basin, and that of a circle with the same perimeter:

$$C_B = 4\pi A/P^2,$$

Basin Compactness is the perimeter to area ratio:

$$C = P/A,$$

The Melton Ruggedness Number [Mel65] is a measure of the relief in a basin, computed as the dimensionless ratio between the
<figure><img src="image_34.png" /><figcaption>Figure 33: Basin metrics computed on drainage basins of Strahler order 4, range of each metric normalized for clarity.</figcaption></figure>

elevation range and the square root of basin area:

$$MRN(\mathbf{p}) = \frac{\Delta h(\mathcal{D})}{A(\mathbf{p})^{1/2}}$$

The hypsometric curve can be used to characterize different drainage basins. This curve measures the relative cumulative area above a given relative elevation threshold. It is commonly renormalized into the [0, 1]² domain. The area below the curve is known as hypsometric integral, and it measures the relative location of mean elevation within the basin elevation range. This measure has been correlated to the erosional maturity of the basin [Str52, BJ21]. A simplified approximate expression for the hypsometric integral $H_i$ of a drainage basin $\mathcal{D}$ is:

$$H_i = \frac{\bar{h}(\mathcal{D}) - \min h(\mathcal{D})}{\Delta h(\mathcal{D})}$$

In some cases, basin descriptors also take into account its streams. For different methods to delineate streams in a basin, please refer to section 7.1. The Drainage Density is the sum of all stream lengths divided by the basin area. The Relief Ratio [Sch56] is defined as the ratio between the basin height range and the length $L$ of its longest stream $\Delta h(\mathcal{D})/L$ and has higher values on steeper streams. The Form factor is defined as:

$$f(\mathbf{p}) = A(\mathbf{p})/L^2,$$

It has lower values for elongated basins. These descriptors correlate with different types of tectonic activity in the basin [RM21]. For example, lower basin compactness, circularity and form factor values indicate higher tectonic activity, similarly to high values of relief ratio. Figure 33 shows some examples.

## 5.4. Soil wetness

Soil water properties are important descriptors for ecosystems, agriculture and soil management.

The Wetness Index [BK79], also called Topographic Wetness Index or Soil Wetness Index, combines drainage area and local slope:

$$w(\mathbf{p}) = \ln \frac{\tilde{A}(\mathbf{p})}{s(\mathbf{p})}$$

in order to quantify potential soil moisture or flow accumulation due to the topography of the terrain. It reaches higher values when the upslope drainage area is large and the slope is flat.

In order to avoid the dependency of $w$ to the flow accumulation algorithm, [MOA09] proposed the Depth-to-water index, which accounts for the distance and accumulated slope towards a water location (Figure 34). Let $\rho \in [0, 1]$ the coefficient of permeability of the terrain. Given an initial position $\mathbf{p}$ and the path $\{\mathbf{p} = \mathbf{q}_0, \mathbf{q}_1, \dots, \mathbf{q}_n = \mathbf{p}_w\}$ of least cumulative slope towards a water location $\mathbf{p}_w$, the index is defined as:

$$d_{\mathcal{R}}(\mathbf{p}) = \sum_{i=1}^{n} s(\mathbf{q}_i) ||\mathbf{q}_i - \mathbf{q}_{i-1}|| \rho(\mathbf{q}_i)$$

Water features like rivers and lakes are assigned a value of 0, large values correspond to less-saturated, drier soils. Similarly, [NCH*11] proposed a normalization of the height field as the elevation offset between each terrain position and its closest water position $\mathbf{p}_w$. The Height Above Nearest Drainage is defined as:

$$h_{\mathcal{R}}(\mathbf{p}) = h(\mathbf{p}) - h(\mathbf{p}_w)$$

Positions in drainage are assigned an elevation of 0. This representation correlates well with different soil environment classes defined from a set of calibrated thresholds.

<figure><img src="image_35.png" /><figcaption>Figure 34: Wetness and depth-to-water indices.</figcaption></figure>

## 5.5. Erosional indices

The Stream Power Equation describes in its most general form the evolution of a landscape height [WT99]:

$$\frac{\partial h}{\partial t} = u - kA^m s^n$$

$h$ is the elevation of the river, $u$ the uplift rate, $k$ an erodibility coefficient, $A$ the drainage area and $m, n$ constants that depend on the terrain. Assuming flow discharge is proportional to catchment area, different compound indices can be derived from the previous equation [MGNP93].

The Stream Power Index (Figure 35) is a relationship that combines drainage area and local slope in a power-law and correlates with the erosive power of the flow through $\mathbf{p}$:

$$SP(\mathbf{p}) = \tilde{A}(\mathbf{p})^m \cdot s(\mathbf{p})^n$$
<figure><img src="image_36.png" /><figcaption>**Figure 35:** Stream power indices with different exponents to model gully formation (left) and sediment transport (right).</figcaption></figure>

The exponents *m* and *n* are usually obtained empirically for a particular basin or erosional effect. For example, *m* = *n* = 1 results in its simplest form with $SP(\mathbf{p}) = \tilde{A}(\mathbf{p}) \cdot s(\mathbf{p})$ and correlates with locations where gullies can develop [MGL91].

Similarly, [MB86] proposes a Sediment Transport Index correlated with soil loss due to overland flow (sheet and rill erosion processes, excluding gully erosion):

$$STI(\mathbf{p}) \propto \tilde{A}(\mathbf{p})^{0.4} (\arctan s(\mathbf{p}))^{1.3}$$

In soil management and agriculture, the long-term average erosion rate is predicted using the Universal Soil Loss Equation [WS78], an empirical model derived from cropland measurements that accounts for rainfall and runoff, soil properties, slope length and steepness, vegetation cover, crop management and support practices. Among these factors, the slope length and steepness (LS) is commonly computed analytically from the height field. Therefore, assuming other conditions are uniform, erosion could be correlated to this term. A recent review of different formulations of this term is given in [KKRO20].

## 5.6. Analysis

In Computer Graphics, fluvial erosion simulations often need to simulate water flowing on the terrain. Solutions derived from [MKM89], which can be implemented efficiently on graphics hardware [MDH07], approximate the flow routing process by distributing the flow proportionally to the elevation difference similarly to [QBCP91]. In contrast, some recent works on tectonic erosion simulations [CBC*16, CCB*18] favour the faster D8 model to compute the stream power erosion. A rationale for this is the reduction of the complexity of the flow routing algorithm to $O(n)$ instead of $O(n \ln n)$ [CBB19]. Finally, [SPF*23] proposed a large-scale terrain generation approach that relies on the multiple flow direction model described in [Hol94]. Works that specifically address river trajectories and shapes in Computer Graphics [PDG*19, PGCG23] use the drainage area computation, with a threshold, to identify pixels that might be part of a river. This is then cleaned and processed to build a river graph.

Flow regularization processes (depression filling or breaching, hybrid methods) have only been utilized recently in Computer Graphics. Most early terrain generation methods, even erosion simulations, produce many local depressions and flat areas, which penalizes the overall hydrological consistency. [CBC*16] introduced

breaching to Computer Graphics to prevent the formation of lakes in a tectonic erosion method. Scott et al. [SD21] compared the impact of depression filling and breaching in example-based terrain synthesis. Moreover, they proposed interleaving depression breaching steps within the generation process to obtain hydrologically consistent terrains at the price of narrow canyons emergence.

Although the topographic wetness index itself has been seldom used in Computer Graphics, many ecosystem simulations and procedural vegetation instantiation methods employ user-defined or computed soil moisture maps derived from the slope and the distance to neighbouring rivers [GLCC17, CGG*17, KGG*20].

Artist feedback (see Section 8) informed us that the drainage area $\tilde{A}$ is a good representative of the flow map, which is not surprising. They also showed preference for maps computed using Multiple Flow Direction routing, explaining that D8 looked too straight and unnatural. They expressed being already familiar with wetness maps. However, Stream Power and Sediment Transport Index maps were perceived as a coarser and slope-based version of a flow map.

## 6. Orometry

Orometry is the field of geomorphology that deals with the measurement of mountains and their properties. While the shape or extent of a mountain is usually treated as a fuzzy entity [FWC04], peaks are perfectly identified as their local maxima.

The most obvious measure for a peak or summit $q_i$ is its elevation $h(q_i)$, but, on its own, it fails to incorporate any notion of relative importance: nearby and slightly lower local maxima have similar elevations, although not considered independent summits.

Some authors have used notions of local relief to identify independent peaks. For example, [DW08] classifies a local maximum $q$ as a peak candidate if the elevation drop in a square neighborhood of a certain size is larger than a given threshold: $\Delta h(\Omega(q)) > T$. However, these definitions suffer from the arbitrariness of the domain and threshold choices. Therefore, richer and more informative metrics have been proposed, as described in [Hel05].

### 6.1. Prominence

It is widely used in mountain engineering circles for compiling summit lists and distinguishing between major and minor peaks [Fry87]. Prominence $\pi(q_i)$ is defined as the shortest vertical distance that must be descended from a peak $q_i$ before being able to reascend towards any other higher peak $q_j$. Let $\Gamma_{ij}$ denote a path between peaks $q_i$ and $q_j$, and $\Gamma_{ij}^*$ be any one of the paths that maximize the minimum elevation along its course (Figure 36), the prominence $\pi(q_i)$ of a peak $q_i$ is then defined as (see Figure 37):

<figure><img src="image_37.png" /><figcaption>**Figure 36:** Saddle paths.</figcaption></figure>

$$\pi(q_i) = h(q_i) - \max(\min(h(\Gamma_{ij}^*))) \quad \forall q_j, h(q_j) > h(q_i)$$

An equivalent definition of prominence is the difference in elevation between the peak $q_i$ and the lowest possible contour encircling
<figure><img src="image_38.png" /><figcaption>Figure 37: Notations for prominence and isolation.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Figure 38: Identified peaks and key saddles with two prominence thresholds. Orange: peaks, green: saddles. Black edges represent the connectivity of the Divide tree (Section 7.3).</figcaption></figure>

$\mathbf{q}_i$ but no other higher peak. It follows that there must always exist a saddle $S$ somewhere on the contour, which is identified as the key-saddle of the peak $\mathbf{q}_i$ and denoted as $K_i$. The prominence is then defined as $\pi(\mathbf{q}_i) = h(\mathbf{q}_i) - h(K_i)$. By definition, the sea level is used as key saddle for island and continental high points.

Note that by inverting the sign of the height field the same algorithm used to compute prominence can output the anti-prominence, which measures how deep a depression/hole is. In this case, the key-saddle is the outlet or outflow point of the hole.

Given two peaks of equivalent prominence, the one with a higher prominence to elevation ratio is perceptually more significant. This ratio is called dominance:

$$d(\mathbf{q}_i) = \pi(\mathbf{q}_i) / h(\mathbf{q}_i)$$

Jurgalski [Jur16] employs this metric to distinguish between peaks, sub-peaks, mountain high points, and massif high points. Figure 38 illustrates the identification of peak and saddles using different prominence thresholds.

## 6.2. Isolation

The isolation of a peak (Figure 37) is the minimum Euclidean distance from $\mathbf{q}_i$ to any higher point on the terrain:

$$i(\mathbf{q}_i) = \min_{\mathbf{p} \in \Omega, h(\mathbf{p}) > h(\mathbf{q}_i)} \|\mathbf{p} - \mathbf{q}_i\|$$

From a viewing perspective, the summits of isolated mountains offer the best vantage points since no other peaks occlude the vista. Only the highest point of a given domain has undefined isolation, e.g. Mount Everest on Earth.

While isolation is given on a per-peak basis, a related metric is the Nearest Neighbor Index of a domain $\Omega$ of planar area

$A_{\Omega}$ [CE54]. Consider the set of $n$ peaks in $\Omega$, let $d(\mathbf{q}_i)$ denote the distance between $\mathbf{q}_i$ and the nearest peak:

$$d(\mathbf{q}_i) = \min_{j \neq i} d(\mathbf{q}_i, \mathbf{q}_j)$$

The Nearest Neighbor Index is defined as the ratio between the mean distance to nearest neighbor and the expected mean distance in a random distribution of locations:

$$NNI(\Omega) = \left( \frac{1}{n} \sum_{\mathbf{q}_i \in \Omega} d(\mathbf{q}_i) \right) / \left( \frac{1}{2} \sqrt{A_{\Omega}/n} \right)$$

Values are close to 0 for clustered distributions, to 1 for randomly distributed peaks, and to 2.15 for a regular distribution forming an equilateral triangle tiling.

## 6.3. Impressiveness

None of the metrics defined so far account for the shape of the mountains. Climbers often prefer impressive peaks, and while this is a highly subjective matter, some efforts have been made at proposing descriptors correlated with this perception.

Omni-directional Relief and Steepness [EM15] accounts for both elevation and steepness of a mountain (Figure 39). When comparing two peaks of similar altitude, the steepest one has a larger value; in the case of two peaks with similar steepness, the highest one has a larger value. Mathematically, it is the weighted average of directional slopes close to a root mean square:

$$ORS(\mathbf{p}) = \left[ \int_{\Omega} f(s(\mathbf{p}, \mathbf{q}))^2 d\mathbf{q} \right]^{1/2}$$

The function $f$ is chosen to be a monotonically increasing function that ensures integral convergence. The authors suggest using:

$$f(u) = \left[ \frac{1}{\pi^3} (2u \arctan u - \ln(u^2 + 1) - \arctan^2 u) \right]^{1/2}$$

Inspired by ORS, the Jut metric (Figure 40) proposed by [Xu22] defines the impressiveness of a feature as the maximum angle-reduced height of $\mathbf{p}$ with respect to any other point $\mathbf{q}$:

$$j(\mathbf{p}) = \max_{\mathbf{q} \in \Omega} z_{\mathbf{p} \leftarrow \mathbf{q}} \cdot \sin \theta_{\mathbf{p} \leftarrow \mathbf{q}}$$

where $z_{\mathbf{p} \leftarrow \mathbf{q}}$ is the height of $\mathbf{p} \in \mathbb{R}^3$ above the horizontal plane (perpendicular to gravity) at point $\mathbf{q} \in \mathbb{R}^3$, and $\theta_{\mathbf{p} \leftarrow \mathbf{q}}$ is the angle of elevation of $\mathbf{p}$ with respect to $\mathbf{q}$. Jut reaches larger values - more impressiveness - when the terrain rises abruptly in a small area, such as in the vicinity of cliffs.

Apart from jut, Xu's work also proposes several metrics related to impressiveness that he calls dominance (the maximum height of $\mathbf{p}$ above the horizontal plane of any other point $\mathbf{q}$), submission (the maximum height from $\mathbf{p}$ of any other point $\mathbf{q}$), rut (the maximum angle-reduced height from $\mathbf{p}$ to any other point $\mathbf{q}$), and other composite measures between them. Moreover, their formulation takes into account planetary curvature, so $z_{\mathbf{p} \leftarrow \mathbf{q}}$ can be negative if $\mathbf{p}$ lies below the horizon plane of $\mathbf{q}$. Please refer to [Xu22] for more definitions and details.
<figure><img src="image_40.png" /><figcaption>**Figure 39:** Omni-directional relief and steepness computed at 1 km cut-off highlights local sharp features, whereas the same metric at 10 km cut-off favors the highest landforms.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>**Figure 40:** Jut highlights steep and high features in a landscape, whereas Rut identifies the locations from where such features might appear impressive.</figcaption></figure>

## 6.4. Peakedness

It is also possible to question whether all local maxima (summits) are to be considered peaks. While the mountaineering community has settled for topographic prominence as an inclusion criterion - for example, 30m are required in the list of 4,000m peaks of the Alps, and 300 feet in the *fourteeners* of the United States - the highest point of a large plateau with almost negligible slope could qualify as a peak according to an appropriate prominence threshold.

<figure><img src="image_42.png" /><figcaption>**Figure 41:** Peakedness metric (right) computed as weighted distance to the top 20% peak prototypes in terms of prototypicality measure (left).</figcaption></figure>

Deng [DW08] proposes a fuzzy peakedness metric to grade summits resemblance to prototypical shapes of peaks (Figure 41). It involves the computation of four metrics in the neighbourhood of every summit: elevation range, mean slope, relative elevation of

the summit in a three times larger neighbourhood, and the number of other peaks in this neighbourhood. The peak prototypicality measure, a weighted combination of those metrics results, is used to locate the top-scoring summits and average their scores into a peak prototype descriptor. The fuzzy peakedness of a summit is then given by weighting its metric distances to those of the peak prototypes using various neighbourhood sizes.

## 6.5. Analysis

Orometry, particularly prominence and isolation, was successfully introduced in Computer Graphics for analyzing and generating terrains in [AGP*19]. By constructing a sparse metric computed from the peaks and saddles, it is possible to robustly distinguishing between different types of mountain ranges. The synthesis algorithm builds a graph of peaks and saddles respecting a given orometric distribution and incrementally derives a consistent river network as a dual graph before shaping the valley.

The village generation algorithm described in [EBP*12] somehow used metrics that resemble prominence to identify locations with favorable conditions to determine where new villages could develop. The proposed algorithm combines slope $s(p)$, accessibility $\alpha(p)$, distance to water bodies and other factors, particularly the so-called geographical domination defined as:

$$f(\mathbf{p}) = \frac{1}{A_{\Omega}} \int_{\Omega} \frac{h(\mathbf{p}) - h(\mathbf{q})}{1 + \|\mathbf{p} - \mathbf{q}\|^2} d\mathbf{q}$$

Geographical domination for a given point $p$ is negative if its elevation is below the average elevation of the considered region $\Omega$, and is the higher as its elevation is greater than the surroundings. Note that Viewsheds (Section 3.1) could be employed and combined with other metrics to identify regions of particular importance for village defense.

## 7. Surface networks

From a geomorphological perspective, a terrain is characterized by two interlocking networks with a dendritic tree structure: the channel network, which represents streams and rivers, and the ridge network, which connects peaks and delineates watersheds [Wer88].

The 19th century works by Cayley [Cay59] and Maxwell [Max70] were the first to describe a terrain surface in terms of peaks, pits, saddles, ridge lines and course lines. In particular, Maxwell detailed a procedure similar to the watershed transform for identifying these features. Starting at the elevation level of the lowest pit, this level value is progressively raised resulting in different contour curve sets. Initially, a closed curve appears around the pit, bounding the depression around it. More pits and their depressions successively appear as the level is raised. When two contours corresponding to two distinct depressions merge into a point, this location is called a bar and the merging depressions are labeled as one (Figure 42 left). In contrast, when two points on the same depression contour merge, two distinct regions of higher elevation than the level set become disconnected, each with its respective local maximum (peak): this point corresponds to a pass (Figure 42 right). There are $p-1$ passes for $p$ peaks and $q-1$ bars for $q$ pits. Note that both bars and passes correspond to
<figure><img src="image_43.png" /><figcaption>Figure 42: Maxwell's procedure to locate bars and passes.</figcaption></figure>

saddle points, i.e. non-extrema critical points. The steepest ascent lines from a saddle to neighboring peaks delineate the ridge lines, and the steepest descent lines the course lines. Moreover, Maxwell defined hills and dales as, respectively, the area around a peak encircled by course lines, and the area around a pit enclosed by ridge lines.

Pfaltz [Pfa76] was the first to formally define a surface network and its topological abstraction as a graph. He extended Morse theory by incorporating Maxwell's ideas in order to build the connectivity between the surface's critical points through the ridge and course lines. Thus, a surface network can be seen as a simplified and discrete representation of the topology of a Morse-Smale complex. We refer the reader to [CR16] for a review of surface network algorithms and their applications, and to Chapter 2 in [RM02] for a discussion on their evolution and alternative graph representations.

## 7.1. Channel networks

The channel network – also called river network, stream network or drainage network – is a graph formed by the set of river edges and their joint nodes (Figure 43). Formerly defined manually, it is common nowadays to extract it automatically from the terrain height field [Lin06] by computing the drainage area (see Section 5.1) and selecting cells above a threshold value [OM84]: $R = \{\mathbf{p}_i | \tilde{A} \cdot s^2 > C\}$.

This approach produces unrealistic drainage networks with uniform density. The Channel Initiation Threshold introduced in [MFG93] addresses this limitation by accounting for the local slope such that the drainage area needed to form a channel is smaller in steeper areas:

$$R = \{\mathbf{p}_i | \tilde{A} \cdot s^2 > C\}$$

The threshold C depends on the surface properties and the climate.

For a given basin, the river network takes the form of a binary tree: in nature, it is exceptional to find confluences of more than two rivers [Shr66]. The root is the outlet or river mouth, the lowest point downstream, and the leaves correspond to the sources, the furthest points upstream. Connections or graph edges are channel sections between a source or a river fork and the next fork downstream. Each node contains two tributary links and one downstream link. Lakes and dams are often ignored and simplified to links.

The classical work by Horton [Hor45] proposed an ordering scheme for channel networks. Strahler [Str57] later simplified it to the following topological rules: 1) all source streams have order

<figure><img src="image_44.png" /><figcaption>Figure 43: Channel network computed from D8 stream area computation. River width is drawn proportionally to drainage area.</figcaption></figure>

1; 2) when two streams of the same order $\omega$ join, the downstream is assigned order $\omega + 1$; 3) when two streams of different order join, the highest order is propagated downstream (Figure 44). Let $\omega_d$ denote the order of the downstream link and $\{\omega_u\}$ the orders of incoming upstream links respectively, Strahler's ordering writes as:

$$\omega_d = \max(\min\{\omega_u\} + 1, \max\{\omega_u\})$$

Horton's ordering included an additional geometric rule to propagate the orders upstream, following the tributary more parallel to the merged stream, or the longest one in case of similar angles. By definition, the order of a basin is the highest order of its streams.

<figure><img src="image_45.png" /><figcaption>Figure 44: Horton-Strahler ordering.</figcaption></figure>

The Horton-Strahler ordering has been used in hydrology to derive properties about channel networks. For example, Horton's law of stream numbers states that

the number $N_\omega$ of streams of order $\omega$ in a network follows an inverse geometric progression with constant bifurcation ratio $R_B$, Horton's law of stream lengths states that average stream length $\tilde{L}_\omega$ increases geometrically with length ratio $R_L$, Hack's law [Hac57] relates to a power law the main stream length and the drainage basin area, and Tokunaga [Tok66] counted the number $T_{n,m}$ of side tributaries of order $m$ into a stream of order $n > m$ to derive additional topological properties. Abrahams [Abr84] provides a review of stochastic channel network models and derived laws.

The Sinuosity of a river is commonly defined as the ratio of its length divided by the linear distance between its end points. Mueller [Mue68] proposed two indices to differentiate between sinuosity caused by the adjustment to terrain topography and sinuosity caused by hydraulic effects, such as meandering. Let $L_s$ denote the stream curve length, $L_v$ the length of a curve placed midway between

<figure><img src="image_46.png" /><figcaption>Figure 45: Sinuosity and channel metrics.</figcaption></figure>
<figure><img src="image_47.png" /><figcaption>Figure 46: Rosgen classification of rivers.</figcaption></figure>

valley walls (or the average between the curve lengths of the left and right valley sides, see Figure 45), and d the distance between the source and mouth of the stream, [Mue68] defines the Channel Index ($c_R$) and Valley Index ($v_R$) as:

$$c_R = \frac{L_s}{d} \qquad VI = \frac{L_v}{d}$$

The Hydraulic Sinuosity Index ($h_R$) and Topographic Sinuosity Index ($t_R$) are derived as:

$$h_R = \frac{c_R - v_R}{c_R - 1} \qquad t_R = \frac{v_R - 1}{c_R - 1} = 1 - h_R$$

Note that those indexes are complementary. In young streams confined by valley walls, topographic sinuosity dominates, whereas, in mature floodplains, hydraulic effects replace those of topography due to lateral migrations and meandering.

## 7.2. River classifications

Stream order can be used as a systematic categorization of the streams in a channel network. Strahler number correlates with magnitudes like stream length or drainage area. However, depending on how the channel network is defined, the number of streams, and thus, the ordering, may vary. Moreover, the Strahler number does not provide information about the geological processes that shaped the rivers, and streams of the same order in different basins can be radically different.

There exist many river classification schemes based on other criteria, qualitative or quantitative, such as the sediment transport type, the geometric pattern of the channel, the river bed material, or the different geomorphic units that can be found in the channel. For a recent review, we refer the reader to [BM22]. One of the most used classifications in land management and stream restoration is Rosgen's classification of natural rivers [Ros94]. He proposed eight major stream types based on the number of channels (multiple or single), the entrenchment ratio (width of the flood-prone area to the river bank width), the channel width-to-depth ratio, and its sinuosity (see Figure 46). These classes are further refined into a total of 94 subclasses depending on slope and bedrock type.

Montgomery and Buffington [MB97] distinguish three kinds of valley substrates: bedrock, colluvial and alluvial. The latter is further subdivided into different morphologies: cascade, step-pool, plane-bed, pool-riffle, and dune-ripple. Although the identification is visual, each class correlates with physical parameters like slope, drainage area, sediment grain size, and channel transport capacity.

Sediment transport shapes the morphology of the river. Schumm [Sch85] classifies alluvial channels into five broad patterns and explains their morphology and characteristics in sediment type, width-to-depth ratio, channel gradient, and stability. Church [Chu06] builds on that classification relating these properties to Shields number, a dimensionless metric of the river shear stress that quantifies its ability to displace sediments of a given size.

## 7.3. Ridge networks

Ridge networks, sometimes called divide networks, act as division lines between drainage basins. Like river networks, they also exhibit a dendritic tree-like structure. Cycles may exist under the presence of endorheic basins, i.e. enclosed depressions that do not flow to a sea or terrain boundary, but are usually removed using a filling/breaching process for hydrological analysis. Contrary to rivers, which are monotonically descending in elevation, such lines have jagged profiles that go up and down with intermediate peaks. Moreover, although confluences of three or more rivers at one location are infrequent, it is typical that several ridge lines converge to one peak. Figure 47 compares several segmentation strategies.

In a similar spirit as done with river networks, it is possible to extract ridge networks by thresholding or applying morphological operations on some attributes such as the Laplacian or curvatures. A low Drainage Area might also be used as an indicator of ridges when computed using a converging algorithm like D8. However, these approaches usually lead to unconnected results and require tuning thresholds or some merging procedure if a continuous network is desired.

The Profile recognition and Polygon-breaking Algorithm (PPA) [CSH98], designed for height fields, often yields better results and proceeds as follows. First, each cell is marked as a ridge candidate if lower neighbors within a distance n can be found on both sides of an axis (most often the cardinal directions ±x, ±y). Then, neighboring candidate points are connected to form a graph with edges weighted by the average elevation of the two endpoints. The original algorithm then visits all segments in increasing elevation and removes them if they belong to a cycle. This computationally demanding approach was improved by using a Minimum-Spanning Tree algorithm instead in [BdBG10]. Eventually, small spurious terminal branches are pruned if their length is smaller than a threshold, typically half the profile length.

Following the ideas on surface network construction [Max70],
<figure><img src="image_48.png" /><figcaption>**Figure 47:** Ridge segmentation comparison: drainage area lower than 1 cell, minimum curvature threshold, PPA, ridges from divide tree.</figcaption></figure>

Earl [Ear16] developed an exact, deterministic, and parameter-free algorithm to detect the ridge network. A modified implementation [KdF17] was employed to compute the prominence of the mountains on Earth from a global Digital Elevation Model. First, the algorithm labels all peaks and saddles. A peak is a cell or continuous flat region whose boundary cells have a lower elevation. A saddle is a cell or a region whose boundary cells contain at least two disjoint higher sections. If a saddle boundary contains $N > 2$ higher sections, it is logically considered as $N - 1$ saddles with exactly $N = 2$. The steepest ascent lines are followed from the two higher regions of each saddle until a peak is reached. Computing the prominence of the peaks requires the computation of the key saddles only and yields the *Divide Tree* $T_D$, which is a subset of the ridge network (see Figure 38 for examples). This representation removes the lowest saddle of each closed loop and its corresponding two ridges. Such saddles represent the outflow points of closed basins and were consequently named *Basin Saddles*, or bars [Max70]. Pits or depressions are rare in the case of terrains shaped by fluvial erosion. One notable limitation of the Divide Tree is that it does not contain decreasing lines from a peak that do not end in saddles. Moreover, it binds the peaks into a single connected graph: the Divide Tree always connects two peaks following the least descent path even if no geomorphological ridge exists between them (e.g., groups of volcanoes or dunes).

## 7.4. Ridge hierarchies

One significant difference between channel and ridge networks arises when assigning a hierarchy or order to ridges. While straightforward for channel networks because streams form directed descending paths towards a downstream junction or outlet, no directionality or monotonicity exists for ridge lines. Moreover, nodes in ridge networks (peaks) can have a valence larger than two: several ridges can meet at the same point.

Rack *et al.* [RKO*19] define the main ridge as the longest path or network diameter (Figure 48). Then, the longest unvisited path starting at each bifurcating node is a second-level ridge, and so on recursively. This approach generates a hierarchy based on path length, which allows for studying the different properties of the network structure of mountain ranges. From a geomorphological perspective however, the longest path might not coincide with the main drainage divide (MDD). For instance, Jancewicz *et al.* [JRS*22] manually delineate the MDD from the computed borders of all the drainage basin boundaries and study several metrics along this MDD, asymmetries between metrics of surface and basins on oppo-

<figure><img src="image_49.png" /><figcaption>**Figure 48:** Ridge levels on a DEM of the Central Pyrenees based on the longest ridge criterion by [RKO*19] and the topological ordering scheme by [SS20].</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>**Figure 49:** Maximum Branch Length [LS13] can be used to measure the importance of different drainage divides.</figcaption></figure>

site sides, and surface descriptors in a narrow band close to MDD. We refer to Table 4 in their paper for a summary of metrics and motivation for their use.

Scherler and Schwanghart [SS20] proposed a topographic ordering of ridge trees similar to rivers (Figure 48). Given an acyclic tree-like divide network, the algorithm starts at the endpoint nodes, those with only one incident ridge. These ridges are assigned order 1 and are removed from the divide network. Neighboring nodes that had two incident ridges become endpoints with only one remaining ridge, which will be assigned order 2 and will be removed. This procedure iterates until all ridges have been ordered. Network directionality can also be obtained by considering the last remaining ridge as the outgoing direction from the node, and all the others as incoming ones. The obtained ordering follows a topographical scheme: the degree $\omega_o$ of the outgoing ridge is one unit higher than
<figure><img src="image_51.png" /><figcaption>**Figure 50:** Examples of terrains covered with rocks, snow and vegetation; vegetation density is computed from abiotic parameters such as slope, direct and indirect lighting, whereas snow and rocks are distributed according to slope and flow maps.</figcaption></figure>

the largest degree of incoming ridges $ω_o = \max(\{ω_i\}) + 1$. Finally, note that order is assigned to segments, whereas [RKO*19] labeled entire ridge paths.

Lindsay and Seibert [LS13] proposed a metric based on the convergence path length between flows at opposite sides of the mountain ridges to characterize the importance of different divide ridges. The method relies on the observation that flow paths starting at opposite sides of the divide probably converge after a small distance for small or secondary ridges. In contrast, the flow paths extend over large distances or may not converge for the main divides. The *Branch length* $B(p, q)$ between two positions **p** and **q** is the distance along the flow path (using the D8 model) initiated at **p** until reaching the confluence with the flow path initiated at **q**. It is not symmetric: $B(p, q) = 0$ if **p** is downstream of **q**. The *Maximum Branch Length* $B_Ω(p)$ computes the longest branch length between **p** and its neighbors **q** ∈ Ω. It can be interpreted as the complement of the contributing area in the same way that the divide network is the complement to the channel network (Figure 49).

## 7.5. Analysis

Given their close relationship with shape analysis and geometry processing, it is not surprising to find methods in the Computer Graphics literature addressing the computation of surface networks (see [TIS*95] for an early example) and, more generally, Morse-Smale complexes (see [DFFIM15] for a review). These approaches have, in turn, contributed to advances in geomorphology and GIS.

Classifications are often successfully used in Computer Graphics as they define a limited set of archetypes, patterns or configurations that lend themselves for procedural generation. A typical example is the *Rosgen classification* (see Section 7.2) that was first employed as a proxy for defining the geometry of river channels according to their classification based on slope and drainage [GGG*13]. By combining power laws correlating the flow, hence the width and depth of the river to the drainage area A, the authors propose a procedural approach for generating the profile of the riverbed. Peytavie *et al.* [PDG*19] introduce a similar algorithm to animate the surface of the water. The *Horton Strahler numbering* (see Section 7.2) is also used for classifying rivers into different categories.

Paris *et al.* [PGCG23] modeled the evolution of meandering rivers, using sinuosity and upstream curvature to compute the migration rate along the channel. For an overview of the different methods available to model rivers, see the work of Brown [BP19].

Ridge networks and divide trees (Section 7.3) were employed for landscape classification and synthesis in [AGP*19]. A notable result is the ability to generate a vast variety of terrains conforming to the identified ridge and divide tree patterns. Some image-based generation algorithms such as [ZSTR07] rely on the *Profile recognition and Polygon-breaking Algorithm* to extract ridge networks.

Terrain authoring tools based on sketches typically rely on the user partly delineating rivers and ridges [GMS09, HGA*10, GPM*22]. Learning-based sketching techniques like [GDG*17, LGP*23] can benefit from robust algorithms and metrics for detecting and filtering these networks.

Network representations have also been used in other procedural modeling works apart from terrains. For instance, AlHalawani *et al.* [AYWM14] extracted topological and geometrical properties of street layouts to predict their city or region, and then proposed a simulated annealing optimization to modify an initial random layout into a desired city style. Polasek *et al.* [PHBv21] also analyze the branching structure of trees to learn a realism metric of tree models.

Artist feedback on the maximum branch length $B_Ω$ tells us it has a potential to detecting mountain ridges. Detecting such features is important to place materials such as rocks and in the same time avoid populating with vegetation that is not supposed to grow in these areas.

## 8. Discussion

Here we conclude this work by providing Computer Graphics artists feedback on the different previously described metrics, we compare their computation cost and correlations, and propose future research directions.

**Artists feedback** To assemble a preliminary and qualitative study of the impact the different metrics could have on production
<table><caption>**Table 1:** Comparison table, part 1. P means pointwise, F featurewise, ○ means not used, ● intensively used and ● is intermediate. Cost factor relative to average slope time computation (18 ms), metrics defined on Ω<sub>n</sub> or having cut-off distance computed for n = 1, 5, 16, 53.</caption><thead><tr><th>Section</th><th>Notation</th><th>Description</th><th>Type</th><th>Domain/Feature</th><th>Cost factor</th><th>CG usage</th></tr></thead><tbody><tr><td rowspan="5">2.1</td><td>∇h(p)</td><td>Gradient</td><td rowspan="5">P</td><td rowspan="5">Ω<sub>1</sub></td><td rowspan="5">1</td><td>●</td></tr><tr><td>s(p)</td><td>Slope</td><td>●</td></tr><tr><td>s(p, q), s̄(p)</td><td>Directional and average slope</td><td>○</td></tr><tr><td>n(p)</td><td>Normal</td><td>●</td></tr><tr><td>α(p)</td><td>Terrain aspect</td><td>○</td></tr><tr><td rowspan="2">2.2</td><td>Δh(p)</td><td>Laplacian</td><td>P</td><td>Ω<sub>1</sub></td><td>1</td><td>●</td></tr><tr><td>(-Δ)<sup>s</sup>h(p)</td><td>Fractional Laplacian</td><td>P</td><td>Ω<sub>n</sub></td><td>3 52 480 830</td><td>○</td></tr><tr><td rowspan="2">2.3</td><td>κ<sub>min</sub>, κ<sub>max</sub>, κ<sub>G</sub></td><td>Min, max, and Gaussian curvatures</td><td>P</td><td>Ω<sub>1</sub> or Ω<sub>n</sub></td><td>4 23 175 1800</td><td>●</td></tr><tr><td>κ<sub>pro</sub>, κ<sub>con</sub>, κ<sub>tan</sub></td><td>Profile, contour, and tangential curv.</td><td>P</td><td>Ω<sub>1</sub> or Ω<sub>n</sub></td><td></td><td>●</td></tr><tr><td>2.4</td><td>z(x, y)</td><td>Fitted quadric surface</td><td>P</td><td>Ω<sub>1</sub> or Ω<sub>n</sub></td><td>3 21 175 1800</td><td>●</td></tr><tr><td rowspan="6">2.5</td><td>Δh<sub>Ω<sub>n</sub></sub>(p)</td><td>Local relief</td><td rowspan="6">P</td><td rowspan="6">Ω<sub>n</sub></td><td>1 5 50 530</td><td>○</td></tr><tr><td>h̄(p), σ<sup>2</sup><sub>z</sub>(p)</td><td>Local mean and variance</td><td>1 5 30 340</td><td>○</td></tr><tr><td>t(p) and t̂(p)</td><td>Topographic Position Index (TPI)</td><td>1 4 27 300</td><td>○</td></tr><tr><td>r(p) and r̄(p)</td><td>Terrain Ruggedness Index (TRI)</td><td>1 4 34 400</td><td>○</td></tr><tr><td>ρ(p)</td><td>Rugosity</td><td>3 9 65 710</td><td>○</td></tr><tr><td>σ<sub>s</sub>(p)</td><td>Roughness</td><td>6 12 78 960</td><td>○</td></tr><tr><td rowspan="2">2.6</td><td>HA(p)</td><td>Hillslope asymmetry</td><td>P</td><td>Ω<sub>n</sub></td><td>12 18 120 950</td><td>○</td></tr><tr><td>Asymmetry(X)</td><td>Asymmetry index</td><td>F</td><td>Line, Curve</td><td>-</td><td>○</td></tr><tr><td rowspan="3">3.1</td><td>V(p)</td><td>Viewshed</td><td>F</td><td>Point</td><td>2 · 10<sup>5</sup></td><td>○</td></tr><tr><td>V(P)</td><td>Cumulative viewshed</td><td>F</td><td>Points</td><td>|P| · 10<sup>5</sup></td><td>○</td></tr><tr><td>v<sub>idx</sub>(p)</td><td>Visibility index</td><td>P</td><td>Ω<sub>∞</sub></td><td>10<sup>11</sup></td><td>●</td></tr><tr><td rowspan="2">3.2</td><td>Φ<sub>d</sub>(p), Ψ<sub>d</sub>(p)</td><td>Zenith and Nadir angle</td><td rowspan="2">P</td><td rowspan="2">Ω<sub>∞</sub></td><td>2 3 5 16</td><td>○</td></tr><tr><td>Φ<sub>H</sub>(p), Ψ<sub>H</sub>(p)</td><td>Positive and negative openness</td><td>8 20 48 130</td><td>○</td></tr><tr><td rowspan="4">3.3</td><td>φ<sub>d</sub>(p)</td><td>Horizon angle</td><td rowspan="4">P</td><td rowspan="4">Ω<sub>∞</sub></td><td>180</td><td>○</td></tr><tr><td>a(p, n)</td><td>Accessibility</td><td>150 · n<sub>rays</sub></td><td>●</td></tr><tr><td>v<sub>H</sub>(p)</td><td>Sky view factor</td><td>120 · n<sub>rays</sub></td><td>●</td></tr><tr><td>v̄<sub>H</sub>(p)</td><td>Sky view factor approximation</td><td>2</td><td>●</td></tr><tr><td rowspan="2">3.4</td><td>E(p)</td><td>Surface irradiance</td><td>P</td><td>Ω<sub>∞</sub></td><td>2800</td><td>●</td></tr><tr><td>θ(p)</td><td>Diurnal Anisotropic Heat Index</td><td>P</td><td>Ω<sub>1</sub></td><td>5</td><td>○</td></tr><tr><td rowspan="2">4.3</td><td>p(p)</td><td>Percentile</td><td>P</td><td>Ω<sub>n</sub></td><td>1 4 37 360</td><td>○</td></tr><tr><td>W<sup>ε</sup>, B<sup>ε</sup></td><td>White and Black Top Hat</td><td>P</td><td>Ω<sub>n</sub></td><td>6 30 170 1850</td><td>○</td></tr><tr><td>4.5</td><td>Δ<sub>d</sub></td><td>Geomorphon pattern</td><td>P</td><td>Ω<sub>n</sub></td><td>6 15 41 130</td><td>●</td></tr></tbody></table>

pipelines, we worked in coordination with professional artists in the entertainment industry with extensive experience on landscape modeling (see Figures 50 and 51). They provided us with synthetic terrain maps and we computed the set of terrain metrics. We selected those metrics which result in maps over Ω - the majority of them -, in order to provide the artists with grayscale images representing each one on a normalized range. For metrics that depend on parameters, we generated several maps with two or three values of each parameter, providing more variety.

The artists carefully analyzed the different maps, and presented us with comments on each metric. They were able to load them, display them on top of the terrain, and even use them in their con-

tent synthesis pipelines. Particularly, they rated the potential of every unexplored metric compared to commonly used ones (such as slope, Laplacian, or accessibility) concerning an original use in a creative landscape generation application. We reported their detailed feedback in the corresponding Analysis subsections, and summarize their main comments below.

Preliminary results demonstrate that some metrics unknown in their existing production pipeline were potentially extremely useful. Artists identified the maximum branch length B<sub>Ω</sub> as interesting to obtain mountain ridges, and some curvature variants κ as good tools to synthesize high-resolution details on bare bedrock elevation maps. They mentioned that openness Φ(p) could have
<table><caption>**Table 2:** Comparison table, part 2. P means pointwise, F featurewise, ○ means not used, ● intensively used and ● is intermediate. Cost factor relative to average slope time computation (18 ms), metrics defined on Ω<sub>n</sub> computed for n = 1, 5, 16, 53. Prominence, dominance and isolation cost factors include computation of the Divide Tree. River metrics report the cost of computing the river network tree, as metric computations on it are negligible.</caption><thead><tr><th>Section</th><th>Notation</th><th>Description</th><th>Type</th><th>Domain/Feature</th><th>Cost factor</th><th>CG usage</th></tr></thead><tbody><tr><td rowspan="7">5.3</td><td>A(p)</td><td>Drainage area</td><td>P</td><td>Ω<sub>∞</sub></td><td>7 (D8) 13 (MFD)</td><td>●</td></tr><tr><td>~A(p)</td><td>Specific Catchment Area</td><td>P</td><td>Ω<sub>∞</sub></td><td>7 (D8) 13 (MFD)</td><td>○</td></tr><tr><td>D</td><td>Drainage basins</td><td>P</td><td>Ω<sub>∞</sub></td><td>7 (D8) + 5</td><td>○</td></tr><tr><td>C<sub>B</sub>, C, f</td><td>Basin circularity, compactness, form factor</td><td>F</td><td>Basin</td><td>-</td><td>○</td></tr><tr><td>MRN(D)</td><td>Melton Ruggedness Number</td><td>F</td><td>Basin</td><td>-</td><td>○</td></tr><tr><td>Δh(D)/L</td><td>Relief ratio</td><td>F</td><td>Basin</td><td>-</td><td>○</td></tr><tr><td>H<sub>i</sub></td><td>Hypsometric Integral</td><td>F</td><td>Basin</td><td>-</td><td>○</td></tr><tr><td rowspan="5">5.4</td><td>w(p)</td><td>Wetness index</td><td rowspan="5">P</td><td rowspan="5">Ω<sub>∞</sub></td><td>9</td><td>○</td></tr><tr><td>SP(p)</td><td>Stream Power Index</td><td>9 (m = n = 1) 12</td><td>●</td></tr><tr><td>STI(p)</td><td>Sediment Transport Index</td><td>11</td><td>○</td></tr><tr><td>d<sub>R</sub>(p)</td><td>Depth-to-water Index</td><td>20</td><td>○</td></tr><tr><td>h<sub>R</sub>(p)</td><td>Height Above Nearest Drainage</td><td>15</td><td>○</td></tr><tr><td rowspan="6">6</td><td>π(q<sub>i</sub>)</td><td>Prominence</td><td>F</td><td>Peak</td><td>9</td><td>●</td></tr><tr><td>d(q<sub>i</sub>)</td><td>Dominance</td><td>F</td><td>Peak</td><td>9</td><td>●</td></tr><tr><td>i(q<sub>i</sub>)</td><td>Isolation</td><td>F</td><td>Peak</td><td>10</td><td>●</td></tr><tr><td>NNI(Ω)</td><td>Nearest Neighbor Index</td><td>F</td><td>Ω<sub>∞</sub></td><td>130</td><td>○</td></tr><tr><td>ORS(p)</td><td>Omni-directional relief and steepness (ORS)</td><td>P</td><td>Ω<sub>n</sub></td><td>3 50 400 3800</td><td>○</td></tr><tr><td>j(p)</td><td>Jut</td><td>P</td><td>Ω<sub>n</sub></td><td>2 22 200 2100</td><td>○</td></tr><tr><td rowspan="5">7.1</td><td>ω<sub>d</sub></td><td>Strahler number</td><td rowspan="5">F</td><td rowspan="5">River</td><td rowspan="5">90</td><td>●</td></tr><tr><td>h<sub>R</sub>, t<sub>R</sub></td><td>Hydraulic and Topographic Sinuosity Index</td><td>●</td></tr><tr><td>R<sub>B</sub>, R<sub>L</sub>, R<sub>A</sub></td><td>Bifurcation, length and drainage ratio</td><td>○</td></tr><tr><td>T<sub>k</sub></td><td>Tributaries number</td><td>○</td></tr><tr><td>PPA</td><td>PPA ridges</td><td>P</td><td>Ω<sub>n</sub></td><td>72 150 315 790</td><td>●</td></tr><tr><td>7.3</td><td>T<sub>D</sub></td><td>Divide Tree</td><td>P</td><td>Ω<sub>∞</sub></td><td>9</td><td>●</td></tr><tr><td rowspan="2">7.4</td><td>B(p, q)</td><td>Branch Length</td><td rowspan="2">P</td><td rowspan="2">Ω<sub>∞</sub></td><td>-</td><td>○</td></tr><tr><td>B<sub>Ω</sub>(p)</td><td>Branch Length max</td><td>75</td><td>○</td></tr></tbody></table>

a promising potential to determine snow patterns. They also reported that the drainage area ~A would be a good flow representative, whereas the fractional Laplacian could be used as a proxy for sediment deposition maps. Eventually, viewsheds V(p) were rated as interesting for delimiting regions, as omnidirectional visibility for placing villages or castles. In contrast, some of the proposed metrics lack the precision to be used in an artistic pipeline: this was the case for hillslope asymmetry HA(p), river distance h<sub>R</sub>(p), or roughness σ<sub>s</sub>(p). A complete evaluation of those metrics was beyond the scope of this paper, and remains an interesting venue of future work.

**Performance comparison** Tables 1 and 2 summarize the relative computational costs of the metrics presented in this paper. Our implementation is accessible at the following link: https://github.com/oargudo/terrain-descriptors.

Execution times were measured using a 1024² heightfield covering an area of 95²km² in the Alps, with timings averaged over five executions. We report a cost factor as the execution time relative to that of computing the slope (18 ms in our machine), as it is the simplest and most commonly used metric in Computer Graphics. Note that this normalization also, in many cases, factors out the influence of terrain grid dimensions; however, global metrics such as visibility or drainage/ridge networks may still depend on terrain size. For metrics computed on a subdomain Ω<sub>n</sub>, we tested values of n corresponding to radii of 100, 500, 1500 and 5000m. As expected, costs exhibit quadratic growth with respect to n.

It is important to note that all metrics were implemented in C++ and run in a single CPU thread, to provide a fairer comparison between those that can easily be parallelized and those that cannot. Our implementations prioritized simplicity and clarity, following the metric definitions closely, since our goal was to provide an insight into their baseline complexity prior to algorithmic optimizations. Since many metrics are computed independently at every position
<figure><img src="image_52.png" /><figcaption>**Figure 51:** Different environments of virtual world: vegetation density is computed from abiotic parameters such as slope, direct and indirect lighting, rocks and snow distributed according to slope and flow maps, road networks generated to conform to slope and curvature.</figcaption></figure>

sition p, a simple parallel dispatch can significantly speed up computation time. In other cases, specific algorithmic developments on a given metric may exist to improve performance. For example, Summed Area Tables instead of a loop can avoid the quadratic cost on *n* whenever metrics require a sum of values over Ω*n*. We tested this approach on the Topographic Position Index and achieved a constant relative cost of 8, thus making it more efficient for *n* > 10.

Finally, Figure 52 shows the correlations between different metrics. In Computer Graphics, especially in procedural modeling and texturing, it is not always necessary to use exact physical quantities. Metric maps are often normalized and rescaled before application. Consequently, it can be interesting to replace a computationally expensive metric with simpler correlated indices. For example, Fractional Laplacian, TPI, or impressiveness metrics like ORS and Jut could serve as proxies for a Visibility Index, at computational costs that are several orders of magnitude smaller.

## 9. Conclusion

In this paper, we tried to provide a comprehensive review of some of the commonly used metrics in geological and Earth sciences to terrain analysis and synthesis in Computer Graphics. Terrain synthesis has been an active area of research for decades, and recent work demonstrate the beneficial impact of research in geomorphology to propose better algorithms.

A research direction worth investigating would consist in anticipating which algorithms in Computer Graphics could be beneficial to geomorphology. Shape analysis, visibility, lighting, rendering and also exaggerated shading represent a huge part of Computer Graphics techniques that would probably help not only in Cartography but also representing terrains.

## Acknowledgments

This work is funded by the project AMPLI ANR-20-CE23-0001 and EOLE ANR-23-CE56-0008, supported by Agence Nationale de la Recherche (France), and is also part of the project SENDA TED2021-129761B-I00, funded by MCIN/AEI/10.13039/501100011033 and the European Union NextGenerationEU/PRTR. O. Argudo also acknowledges a Maria Zambrano fellowship by Ministerio de Universidades (Spain) and the support of the Department of Research and Universities of the Government of Catalonia (2021 SGR 01035). We thank the reviewers for their suggestions to improve this work. Images in Figures 1, 50, and 51 are courtesy of Benoît Martinez, UBISOFT.

## References

[AAC*17] ARGUDO O., ANDUJAR C., CHICA A., GUÉRIN E., DIGNE J., PEYTAVIE A., GALIN E.: Coherent multi-layer landscape synthesis. *The Visual Computer* 33, 6-8 (2017), 1005–1015. 11

[Abr84] ABRAHAMS A. D.: Channel networks: A geomorphological perspective. *Water Resources Research* 20, 2 (1984), 161–188. 18

[AGP*19] ARGUDO O., GALIN E., PEYTAVIE A., PARIS A., GAIN J., GUÉRIN E.: Orometry-based terrain analysis and synthesis. *ACM Transactions on Graphics* 38, 6 (2019). 17, 21

[AJL*21] ACHARYYA A., JALLU R. K., LÖFFLER M., MEIJER G. G., SAUMELL M., SILVEIRA R. I., STAALS F.: Terrain Prickliness: Theoretical Grounds for High Complexity Viewsheds. In *11th International Conference on Geographic Information Science (GIScience 2021) - Part II* (Dagstuhl, Germany, 2021), Janowicz K., Verstegen J. A., (Eds.), vol. 208 of *Leibniz International Proceedings in Informatics (LIPIcs)*, Schloss Dagstuhl – Leibniz-Zentrum für Informatik, pp. 10:1–10:16. 6

[AYWM14] ALHALAWANI S., YANG Y.-L., WONKA P., MITRA N. J.: What makes London work like London? *Computer Graphics Forum* 33, 5 (2014), 157–165. 21

[BA05] BENES B., ARRIAGA X.: Table Mountains by Virtual Erosion. In
<figure><img src="image_53.png" /><figcaption>**Figure 52:** Correlation matrix for several metrics computed on the example terrain. Metrics have been visually grouped by slope-based (Section 2.1), Laplacians and curvatures (Section 2.2, 2.3), relief and roughness (Section 2.5), visibility (Section 3), hydrology (Section 5) and orometry (Section 6).</figcaption></figure>

*Eurographics Workshop on Natural Phenomena* (2005), Poulin P., Galin E., (Eds.), The Eurographics Association, pp. 33–39. 11

[BA09] BÖHNER J., ANTONIĆ O.: Land-surface parameters specific to topo-climatology. In *Geomorphometry*, Hengl T., Reuter H. I., (Eds.), vol. 33 of *Developments in Soil Science*. Elsevier, 2009, pp. 195–226. 7, 8

[BdBG10] BANGAY S., DE BRUYN D., GLASS K.: Minimum spanning trees for valley and ridge characterization in digital elevation maps. In *Proceedings of the 7th International Conference on Computer Graphics, Virtual Reality, Visualisation and Interaction in Africa* (New York, NY, USA, 2010), Association for Computing Machinery, pp. 73–82. 19

[BJ21] BAJRACHARYA P., JAIN S.: Characterization of drainage basin hypsometry: A generalized approach. *Geomorphology* 381 (2021), 107645. 14

[BK79] BEVEN K. J., KIRKBY M. J.: A physically based, variable contributing area model of basin hydrology. *Hydrological Sciences Bulletin* 24, 1 (1979), 43–69. 14

[BLM14] BARNES R., LEHMAN C., MULLA D.: Priority-flood: An optimal depression-filling and watershed-labeling algorithm for digital elevation models. *Computers & Geosciences* 62 (2014), 117–127. 13

[BM22] BUFFINGTON J. M., MONTGOMERY D. R.: Geomorphic classification of rivers: An updated review. In *Treatise on Geomorphology (Second Edition)*, Shroder J. J. F., (Ed.), second edition ed. Academic Press, Oxford, 2022, pp. 1143–1190. 19

[BP19] BROWN R. A., PASTERNACK G. B.: How to build a digital river. *Earth-Science Reviews* 194 (2019), 283–305. 21

[Bro14] BROWN L.: Texture shading: A new technique for depicting ter-
rain relief. In Proceedings of the 9th ICA Mountain Cartography Workshop (2014). 6

[Cay59] CAYLEY A.: On contour and slope lines. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 18, 120 (1859), 264–268. 17

[CBB19] CORDONNIER G., BOVY B., BRAUN J.: A versatile, linear complexity algorithm for flow routing in topographies with depressions. *Earth Surface Dynamics* 7, 2 (2019), 549–562. 13, 15

[CBC*16] CORDONNIER G., BRAUN J., CANI M.-P., BENES B., GALIN E., PEYTAVIE A., ÉRIC GUÉRIN: Large scale terrain generation from tectonic uplift and fluvial erosion. *Computer Graphics Forum* 35, 2 (2016), 165–175. 5, 6, 15

[CCB*18] CORDONNIER G., CANI M.-P., BENES B., BRAUN J., GALIN E.: Sculpting mountains: Interactive terrain modeling based on subsurface geology. *IEEE Transactions on Visualization and Computer Graphics* 24, 5 (2018), 1756–1769. 5, 15

[CE54] CLARK P. J., EVANS F. C.: Distance to nearest neighbor as a measure of spatial relationships in populations. *Ecology* 35, 4 (1954), 445–453. 16

[CEG*18] CORDONNIER G., ECORMIER P., GALIN E., GAIN J., BENES B., CANI M.-P.: Interactive Generation of Time-evolving, Snow-Covered Landscapes with Avalanches. *Computer Graphics Forum* 37, 2 (2018), 497–509. 5, 6, 8

[CGG*17] CORDONNIER G., GALIN E., GAIN J., BENES B., GUÉRIN E., PEYTAVIE A., CANI M.-P.: Authoring Landscapes by Combining Ecosystem and Terrain Erosion Simulation. *ACM Transactions on Graphics* 36, 4 (2017), 134. 8, 15

[Chu06] CHURCH M.: Bed material transport and the morphology of alluvial river channels. *Annual Review of Earth and Planetary Sciences* 34, 1 (2006), 325–354. 19

[CR16] CLARKE K. C., ROMERO B. E.: On the topology of topography: a review. *Cartography and Geographic Information Science* 44, 3 (2016), 271–282. 18

[CSH98] CHANG Y.-C., SONG G.-S., HSU S.-K.: Automatic extraction of ridge and valley axes using the profile recognition and polygon-breaking algorithm. *Computers & Geosciences* 24, 1 (1998), 83–93. 19

[CTV*17] CERVILLA A. R., TABIK S., VÍAS J., MÉRIDA M., ROMERO L. F.: Total 3d-viewshed map: Quantifying the visible volume in digital elevation models. *Transactions in GIS* 21, 3 (2017), 591–607. 6

[DF90] DOZIER J., FREW J.: Rapid calculation of terrain parameters for radiation modeling from digital elevation data. *IEEE Transactions on Geoscience and Remote Sensing* 28, 5 (1990), 963–969. 7, 8

[DFFIM15] DE FLORIANI L., FUGACCI U., IURICICH F., MAGILLO P.: Morse complexes for shape segmentation and homological analysis: discrete models and algorithms. *Computer Graphics Forum* 34, 2 (2015), 761–785. 21

[Dik89] DIKAU R.: The application of a digital relief model to land-form analysis in geomorphology. *Three dimensional applications in geographical information systems* (1989), 51–77. 9

[DL22] DAOD M., LAAMRI E. H.: Fractional laplacians: A short survey. *Discrete and Continuous Dynamical Systems – S* 15, 1 (2022), 95–116. 3

[DR95] DUBAYAH R., RICH P. M.: Topographic solar radiation models for gis. *International Journal of Geographical Information Systems* 9, 4 (1995), 405–419. 7

[DW08] DENG Y., WILSON J. P.: Multi-scale and multi-criteria mapping of mountain peaks as fuzzy entities. *International Journal of Geographical Information Science* 22, 2 (2008), 205–218. 15, 17

[Ear16] EARL E.: Winprom. 2016. URL: https://github.com/edwardearl/winprom. 20

[EBP*12] EMILIEN A., BERNHARDT A., PEYTAVIE A., CANI M.-P., GALIN E.: Procedural generation of villages on arbitrary terrains. The

Visual Computer (Proceedings of Computer Graphics International) 26, 6-8 (2012), 809–818. 6, 17

[EM15] EARL E., METZLER D.: Cloud-Capped Towers: Capturing Terrain Characteristics Using Topographic Functionals. *Quaestiones Geographicae* 34, 4 (2015), 7–23. 16

[EPCV15] EMILIEN A., POULIN P., CANI M.-P., VIMONT U.: Interactive procedural modelling of coherent waterfall scenes. *Computer Graphics Forum* 34, 6 (2015), 22–35. 11

[Eva80] EVANS I. S.: An integrated system of terrain analysis and slope mapping. *Zeitschrift fur Geomorphologie* 36 (1980), 274–295. 4

[FB07] FOLDES D., BENES B.: Occlusion-Based Snow Accumulation Simulation. In Workshop in Virtual Reality Interactions and Physical Simulation (2007), Dingliana J., Ganovelli F., (Eds.), The Eurographics Association. 8

[Fis92] FISHER P. F.: First experiments in viewshed uncertainty: simulating fuzzy viewsheds. *Photogrammetric engineering and remote sensing* 58 (1992), 345–345. 6

[FL91] FAIRFIELD J., LEYMARIE P.: Drainage networks from grid digital elevation models. *Water Resources Research* 27, 5 (1991), 709–717. 12

[Flo16] FLORINSKY I.: *Digital terrain analysis in soil science and geology*. Academic Press, 2016. 4

[Fry87] FRY S.: Defining and sizing-up mountains. *Summit, Jan.-Feb.* (1987), 16–21, 32. 15

[FSH20] FARMAKIS-SEREBRYAKOVA M., HURNI L.: Comparison of relief shading techniques applied to landforms. *ISPRS International Journal of Geo-Information* 9, 4 (2020). 6, 8

[FW98] FISHER P., WOOD J.: What is a mountain? Or the englishman who went up a boolean geographical concept but realised it was fuzzy. *Geography* (1998), 247–256. 10

[FWC04] FISHER P., WOOD J., CHENG T.: Where is helvellyn? Fuzziness of multi-scale landscape morphometry. *Transactions of the Institute of British Geographers* 29, 1 (2004), 106–128. 10, 15

[GDG*17] GUÉRIN E., DIGNE J., GALIN E., PEYTAVIE A., WOLF C., BENES B., MARTINEZ B.: Interactive example-based terrain authoring with conditional generative adversarial networks. *ACM Transactions on Graphics (Proceedings of Siggraph Asia 2017)* 36, 6 (2017). 1, 21

[GDGP16] GUÉRIN E., DIGNE J., GALIN E., PEYTAVIE A.: Sparse representation of terrains for procedural modeling. *Computer Graphics Forum* 35, 2 (2016), 177–187. 11

[GGG*13] GÉNEVAUX J.-D., GALIN É., GUÉRIN É., PEYTAVIE A., BENEŠ. B.: Terrain generation using procedural models based on hydrology. *ACM Transaction on Graphics* 32, 4 (2013), 143:1–143:13. 11, 21

[GGP*15] GÉNEVAUX J.-D., GALIN E., PEYTAVIE A., GUÉRIN E., BRIQUET C., GROSBELLET F., BENES B.: Terrain Modelling from Feature Primitives. *Computer Graphics Forum* 34, 6 (2015), 198–210. 11

[GGP*19] GALIN E., GUÉRIN E., PEYTAVIE A., CORDONNIER G., CANI M.-P., BENES B., GAIN J.: A review of digital terrain modeling. *Computer Graphics Forum* 38, 2 (2019), 553–577. 1, 5

[GLCC17] GAIN J., LONG H., CORDONNIER G., CANI M.-P.: Eco-brush: Interactive control of visually consistent large-scale ecosystems. *Computer Graphics Forum* 36, 2 (2017), 63–73. 8, 15

[GM97] GARBRECHT J., MARTZ L. W.: The assignment of drainage direction over flat surfaces in raster digital elevation models. *Journal of Hydrology* 193, 1 (1997), 204–213. 13

[GMS09] GAIN J., MARAIS P., STRASSER W.: Terrain sketching. In *Proceedings of the 2009 Symposium on Interactive 3D Graphics and Games* (New York, NY, USA, 2009), I3D ’09, Association for Computing Machinery, p. 31–38. 21
[GPM*22] GUERIN E., PEYTAVIE A., MASNOU S., DIGNE J., SAUVAGE B., GAIN J., GALIN E.: Gradient Terrain Authoring. Computer Graphics Forum 44, 2 (2022), 85–95. 6, 21

[GPMG10] GALIN E., PEYTAVIE A., MARÉCHAL N., GUÉRIN E.: Procedural generation of roads. Computer Graphics Forum 29, 2 (2010), 429–438. 6

[GW00] GALLANT J. C., WILSON J. P.: Primary topographic attributes. Terrain analysis: Principles and applications (2000), 51–85. 3, 4, 10

[Hac57] HACK J. T.: Studies of longitudinal stream profiles in Virginia and Maryland. Tech. rep., 1957. Report. 18

[Hel05] HELMAN A.: The Finest Peaks. Prominence and Other Mountain Measures. Trafford Publishing, 2005. 15

[HGA*10] HNAIDI H., GUÉRIN E., AKKOCHE S., PEYTAVIE A., GALIN E.: Feature based terrain generation using diffusion equation. Computer Graphics Forum 29, 7 (2010), 2179–2186. 21

[Hol94] HOLMGREN P.: Multiple flow direction algorithms for runoff modelling in grid based elevation models: An empirical evaluation. Hydrological Processes 8, 4 (1994), 327–334. 12, 15

[Hor45] HORTON R. E.: Erosional Development of Streams and their Drainage Basins; Hydrophysical Approach to Quantitative Morphology. GSA Bulletin 56, 3 (1945), 275–370. 18

[Jen04] JENNESS J. S.: Calculating landscape surface area from digital elevation models. Wildlife Society Bulletin 32, 3 (2004), 829–839. 4

[JHS*21] JENNY B., HEITZLER M., SINGH D., FARMAKIS-SEREBYAKOVA M., LIU J. C., HURNI L.: Cartographic relief shading with neural networks. IEEE Transactions on Visualization and Computer Graphics 27, 2 (2021), 1225–1235. 6

[JKG*24] JAIN A., KERBL B., GAIN J., FINLEY B., CORDONNIER G.: FastFlow: GPU Acceleration of Flow and Depression Routing for Landscape Simulation. Computer Graphics Forum 43, 7 (2024). 13

[JNS14] JASIEWICZ J., NETZEL P., STEPINSKI T. F.: Landscape similarity, retrieval, and machine mapping of physiographic units. Geomorphology 221 (2014), 104–112. 11

[JRS*22] JANCEWICZ K., RÓZYCKA M., SZYMANOWSKI M., KRYZA M., MIGON P.: Topographic characteristics of drainage divides at the mountain-range scale - A review of dtm-based analytical tools. ISPRS International Journal of Geo-Information 11, 2 (2022), 116. 5, 20

[JS13] JASIEWICZ J., STEPINSKI T. F.: Geomorphons—a pattern recognition approach to classification and mapping of landforms. Geomorphology 182 (2013), 147–156. 10

[Jur16] JURGALSKI E.: Das erweiterte prominenzkonzept. ein mathematisches einteilungssystem für alle berge und gebirge, weltweit anwendbar vom hochgebirge bis zu heimathügeln. Mitteilungen der Fränkischen Geographischen Gesellschaft 61/62 (2016), 105–110. 16

[KdF17] KIRMSE A., DE FERRANTI J.: Calculating the prominence and isolation of every mountain in the world. Progress in Physical Geophysics: Earth and Environment 41, 6 (2017), 788–802. 20

[KGG*20] KAPP K., GAIN J., GUÉRIN E., GALIN E., PEYTAVIE A.: Data-driven Authoring of Large-scale Ecosystems. ACM Transactions on Graphics 39, 6 (2020), 1–14. 6, 8, 15

[KKRO20] KRUK E., KLAPA P., RYCZEK M., OSTROWSKI K.: Influence of dem elaboration methods on the usle model topographical factor parameter on steep slopes. Remote Sensing 12, 21 (2020). 15

[KMN88] KELLEY A. D., MALIN M. C., NIELSON G. M.: Terrain simulation using a model of stream erosion. Computer Graphics 22, 4 (1988), 263–268. 1

[KPJ*21] KENNELLY P. J., PATTERSON T., JENNY B., HUFFMAN D. P., MARSTON B. E., BELL S., TAIT A. M.: Elevation models for reproducible evaluation of terrain representation. Cartography and Geographic Information Science 48, 1 (2021), 63–77. 8

[KS14] KENNELLY P. J., STEWART A. J.: General sky models for illuminating terrains. International Journal of Geographical Information Science 28, 2 (2014), 383–406. 6, 8

[Lan02] LANDIS H.: Production-ready global illumination. 8

[Lea92] LEA N. J.: An aspect driven kinematic routing algorithm. Overland flow: hydraulics and erosion mechanics 147 (1992), 175. 12

[LGP*23] LOCHNER J., GAIN J., PERCHE S., PEYTAVIE A., GALIN E., GUÉRIN E.: Interactive authoring of terrain using diffusion models. Computer Graphics Forum 42, 7 (2023), e14941. 21

[Lin06] LINDSAY J. B.: Sensitivity of channel mapping techniques to uncertainty in digital elevation data. International Journal of Geographical Information Science 20, 6 (2006), 669–692. 18

[Lin16] LINDSAY J. B.: Efficient hybrid breaching-filling sink removal methods for flow path enforcement in digital elevation models. Hydrological Processes 30, 6 (2016), 846–857. 13

[Llo01] LLOBERA M.: Building past landscape perception withgis: Understanding topographic prominence. Journal of Archaeological Science 28, 9 (2001), 1005–1014. 10

[Llo03] LLOBERA M.: Extendinggis-based visual analysis: the concept of visualscapes. International Journal of Geographical Information Science 17, 1 (2003), 25–48. 6

[LNF19] LINDSAY J. B., NEWMAN D. R., FRANCIONI A.: Scale-optimized surface roughness for topographic analysis. Geosciences 9, 7 (2019), 322. 5

[LS13] LINDSAY J. B., SEIBERT J.: Measuring the significance of a divide to local drainage patterns. International Journal of Geographical Information Science 27, 7 (2013), 1453–1468. 20, 21

[Max70] MAXWELL J. C.: On hills and dales. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 40, 269 (1870), 421–427. 17, 19, 20

[MB86] MOORE I. D., BURCH G. J.: Physical basis of the length-slope factor in the universal soil loss equation. Soil Science Society of America Journal 50, 5 (1986), 1294–1298. 15

[MB97] MONTGOMERY D. R., BUFFINGTON J. M.: Channel-reach morphology in mountain drainage basins. GSA Bulletin 109, 5 (1997), 596–611. 19

[MdCPBS14] MICHELON DE CARLI D., POZZER C. T., BEVILACQUA F., SCHETINGER V.: Procedural generation of 3d canyons. In 2014 27th SIBGRAPI Conference on Graphics, Patterns and Images (2014), pp. 103–110. 11

[MDH07] MEI X., DECAUDIN P., HU B.-G.: Fast hydraulic erosion simulation and visualization on GPU. In 15th Pacific Conference on Computer Graphics and Applications (2007), pp. 47–56. 15

[Mel65] MELTON M. A.: The geomorphic and paleoclimatic significance of alluvial deposits in southern arizona. The Journal of Geology 73, 1 (1965), 1–38. 13

[MFG93] MONTGOMERY D. R., FOUFOULA-GEORGIOU E.: Channel network source representation using digital elevation models. Water Resources Research 29, 12 (1993), 3925–3934. 18

[MFS09] MÉNDEZ-FELIU A., SBERT M.: From obscurances to ambient occlusion: A survey. The Visual Computer 25, 2 (2009), 181–196. 8

[MGL91] MOORE I. D., GRAYSON R. B., LADSON A. R.: Digital terrain modelling: A review of hydrological, geomorphological, and biological applications. Hydrological Processes 5, 1 (1991), 3–30. 15

[MGNP93] MOORE I. D., GESSLER P. E., NIELSEN G. A., PETERSON G. A.: Soil attribute prediction using terrain analysis. Soil Science Society of America Journal 57, 2 (1993), 443–452. 14

[MHS*19] MAKOWSKI M., HÄDRICH T., SCHEFFCZYK J., MICHELS D. L., PIRK S., PALUBICKI W.: Synthetic silviculture: Multi-scale modeling of plant ecosystems. ACM Transactions on Graphics 38, 4 (2019). 6, 8

[Mil94] MILLER G.: Efficient algorithms for local and global accessibility shading. In SIGGRAPH (1994), pp. 319–326. 8

[MK02] MCCUNE B., KEON D.: Equations for potential annual direct incident radiation and heat load. Journal of vegetation science 13, 4 (2002), 603–606. 8
[MKM89] MUSGRAVE F. K., KOLB C. E., MACE R. S.: The synthesis and rendering of eroded fractal terrains. *Computer Graphics* 23, 3 (1989), 41–50. 1, 5, 15

[MOA09] MURPHY P. N. C., OGILVIE J., ARP P.: Topographic modelling of soil moisture conditions: a comparison and verification of two models. *European Journal of Soil Science* 60, 1 (2009), 94–109. 14

[MS18] MOKARRAM M., SATHYAMOORTHY D.: A review of landform classification methods. *Spatial Information Research* 26, 6 (2018), 647–660. 9

[Mue68] MUELLER J. E.: An introduction to the hydraulic and topographic sinuosity indexes. *Annals of the Association of American Geographers* 58, 2 (1968), 371–385. 18, 19

[MVNT21] MESTRES N., VERGNE R., NOÛS C., THOLLOT J.: Local light alignment for multi-scale shape depiction. *Computers Graphics Forum* 40, 2 (2021), 575–584. 8

[NCH*11] NOBRE A., CUARTAS L., HODNETT M., RENNÓ C., RODRIGUES G., SILVEIRA A., WATERLOO M., SALESKA S.: Height above the nearest drainage – a hydrologically relevant new terrain model. *Journal of Hydrology* 404, 1 (2011), 13–29. 14

[NWD05] NEIDHOLD B., WACKER M., DEUSSEN O.: Interactive physically based fluid and erosion simulation. In *Proceedings of the Eurographics Workshop on Natural Phenomena* (Dublin, Ireland, 2005), Eurographics Association, pp. 25–32. 1

[OM84] O’CALLAGHAN J. F., MARK D. M.: The extraction of drainage networks from digital elevation data. *Computer Vision, Graphics, and Image Processing* 28, 3 (1984), 323–344. 12, 18

[PD75] PEUCKER T. K., DOUGLAS D. H.: Detection of surface-specific points by local parallel processing of discrete terrain elevation data. *Computer Graphics and image processing* 4, 4 (1975), 375–387. 9

[PDG*19] PEYTAVIE A., DUPONT T., GUÉRIN E., CORTIAL Y., BENES B., GAIN J., GALIN E.: Procedural Riverscapes. *Computer Graphics Forum* 38, 7 (2019), 35–46. 11, 15, 21

[Pfa76] PFALTZ J. L.: Surface networks. *Geographical Analysis* 8, 1 (1976), 77–93. 18

[PGCG23] PARIS A., GUÉRIN E., COLLON P., GALIN E.: Authoring and simulating meandering rivers. *ACM Transactions on Graphics* 42, 6 (2023), 15, 21

[PGP*19] PARIS A., GALIN E., PEYTAVIE A., GUÉRIN E., GAIN J.: Terrain Amplification with Implicit 3D Features. *ACM Transactions on Graphics* 38, 5 (2019), 1–15. 2

[PHBv21] POLASEK T., HRUSA D., BENES B., ČADÍK M.: Ictree: Automatic perceptual metrics for tree models. *ACM Transactions on Graphics* 40, 6 (2021), 21

[PPFB12] POULOS M. J., PIERCE J. L., FLORES A. N., BENNER S. G.: Hillslope asymmetry maps reveal widespread, multi-scale organization. *Geophysical Research Letters* 39, 6 (2012), 10

[PPG*20] PARIS A., PEYTAVIE A., GUÉRIN E., ARGUDO O., GALIN E.: Desertscape Simulation. *Computer Graphics Forum* 38, 7 (2020), 47–55. 6, 11

[PTS99] PREMOŽE S., THOMPSON W. B., SHIRLEY P.: Geospecific rendering of alpine terrain. In *Proceedings of the 10th Eurographics Conference on Rendering* (Goslar, DEU, 1999), EGWR, Eurographics Association, pp. 107–118. 8

[QBCP91] QUINN P., BEVEN K., CHEVALLIER P., PLANCHON O.: The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models. *Hydrological Processes* 5, 1 (1991), 59–79. 12, 15

[RBD06] RUSINKIEWICZ S., BURNS M., DECARLO D.: Exaggerated shading for depicting shape and detail. In *ACM SIGGRAPH 2006* (New York, NY, USA, 2006), SIGGRAPH ’06, Association for Computing Machinery, pp. 1199–1205. 6, 8

[RDE99] RILEY S. J., DEGLORIA S. D., ELLIOT R.: A terrain ruggedness index that quantifies topographic heterogeneity. *Intermountain Journal of sciences* 5, 1–4 (1999), 23–27. 4

[RKO*19] RAK R., KWAPIEŃ J., OŚWIECIMKA P., ZIEBA P., DROŻDŹ S.: Universal features of mountain ridge networks on Earth. *Journal of Complex Networks* 8, 1 (2019), 20, 21

[RKv*22] RAJASEKARAN S. D., KANG H., ČADÍK M., GALIN E., GUÉRIN E., PEYTAVIE A., SLAVÍK P., BENES B.: Ptm: Perceived terrain realism metric. *ACM Transactions on Applied Perception* 19, 2 (2022), 11

[RLD15] REYNOLDS D. T., LAYCOCK S. D., DAY A. M.: Real-time accumulation of occlusion-based snow. *The Visual Computer* 31, 5 (2015), 689–700. 8

[RM02] RANA S. S., MORLEY J.: *Surface Networks*. UCL Working Paper 43, Centre for Advanced Spatial Analysis, University College London, 2002. 18

[RM21] RÓZYCKA M., MIGOŃ P.: Morphometric properties of river basins as indicators of relative tectonic activity – problems of data handling and interpretation. *Geomorphology* 389 (2021), 107807. 14

[RMCRD02] RODRIGUEZ F., MAIRE E., COURJAULT-RADÉ P., DARROZES J.: The black top hat function applied to a DEM: A tool to estimate recent incision in a mountainous watershed (estibère watershed, central pyrenees). *Geophysical Research Letters* 29, 6 (2002), 9–1–9. 4. 10

[Ros94] ROSGEN D. L.: A classification of natural rivers. *CATENA* 22, 3 (1994), 169–199. 19

[RPP93] ROUDIER P., PEROCHE B., PERRIN M.: Landscapes synthesis achieved through erosion and deposition process simulation. *Computer Graphics Forum* 12, 3 (1993), 375–383. 1, 5

[SA05] SCHMIDT J., ANDREW R.: Multi-scale landform characterization. *Area* 37, 3 (2005), 341–350. 10

[SAG98] S. Z., A. I., G. K.: An ambient light illumination model. In *Rendering Techniques* (Wien, 1998), Springer-Verlag, pp. 45–56. 8

[Sch56] SCHUMM S. A.: Evolution of Drainage Systems and Slopes in Badlands at Perth Amboy, New Jersey. *GSA Bulletin* 67, 5 (1956), 597–646. 14

[Sch85] SCHUMM S. A.: Patterns of alluvial rivers. *Annual Review of Earth and Planetary Sciences* 13, 1 (1985), 5–27. 19

[SD21] SCOTT J. J., DODGSON, N. A.: Example-based terrain synthesis with pit removal. *Computers and Graphics* 99, C (2021), 43–53. 15

[SEB03] SCHMIDT J., EVANS I. S., BRINKMANN J.: Comparison of polynomial models for land surface curvature calculation. *International Journal of Geographical Information Science* 17, 8 (2003), 797–814. 4

[SH04] SCHMIDT J., HEWITT A.: Fuzzy land element classification from dtms based on geometry and terrain position. *Geoderma* 121, 3–4 (2004), 243–256. 10, 11

[Sha95] SHARY P. A.: Land surface in gravity points classification by a complete system of curvatures. *Mathematical Geology* 27, 3 (1995), 373–390. 3, 9

[Shr66] SHREVE R. L.: Statistical law of stream numbers. *The Journal of Geology* 74, 1 (1966), 17–37. 18

[Sm14] SMITH M. W.: Roughness in the earth sciences. *Earth-Science Reviews* 136 (2014), 202–225. 5

[SN08] SNYDER J., NOWROUZEZAHRAI D.: Fast soft self-shadowing on dynamic height fields. In *Proceedings of the Nineteenth Eurographics Conference on Rendering* (Goslar, DEU, 2008), EGSR ’08, Eurographics Association, pp. 1275–1283. 8

[Son17] SONI S.: Assessment of morphometric characteristics of chakrar watershed in madhya pradesh india using geospatial technique. *Applied Water Science* 7, 5 (Sep 2017), 2089–2102. 13

[SPF*23] SCHOTT H., PARIS A., FOURNIER L., GUÉRIN E., GALIN E.: Large-scale terrain authoring through interactive erosion simulation. *ACM Transactions on Graphics* 42, 5 (2023), 1–15. 6, 15
[SS20] SCHERLER D., SCHWANGHART W.: Drainage divide networks – part 1: Identification and ordering in digital elevation models. *Earth Surface Dynamics* 8, 2 (2020), 245–259. 5, 20

[SSM02] SHARY P. A., SHARAYA L. S., MITUSOV A. V.: Fundamental quantitative methods of land surface analysis. *Geoderma* 107, 1-2 (2002), 1–32. 4

[SSM05] SHARY P. A., SHARAYA L. S., MITUSOV A. V.: The problem of scale-specific and scale-free approaches in geomorphometry. *Geografia Física e Dinamica Quaternaria* 28, 1 (2005), 81–101. 9

[Str52] STRAHLER A. N.: Hypsometric (Area-Altitude) Analysis of Erosional Topography. *GSA Bulletin* 63, 11 (1952), 1117–1142. 14

[Str57] STRAHLER A. N.: Quantitative analysis of watershed geomorphology. *Eos, Transactions American Geophysical Union* 38, 6 (1957), 913–920. 18

[SVC03] SOILLE P., VOGT J., COLOMBO R.: Carving and adaptive drainage enforcement of grid digital elevation models. *Water Resources Research* 39, 12 (2003). 13

[Tar97] TARBOTON D. G.: A new method for the determination of flow directions and upslope areas in grid digital elevation models. *Water resources research* 33, 2 (1997), 309–319. 12

[TCZR15] TABIK S., CERVILLA A. R., ZAPATA E., ROMERO L. F.: Efficient data structure and highly scalable algorithm for total-viewshed computation. *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing* 8, 1 (2015), 304–310. 6

[TIS*95] TAKAHASHI S., IKEDA T., SHINAGAWA Y., KUNII T. L., UEDA M.: Algorithms for extracting correct critical points and constructing topological graphs from discrete geographical elevation data. *Computer Graphics Forum* 14, 3 (1995), 181–192. 21

[Tok66] TOKUNAGA E.: The composition of drainage network in toyohira river basin and valuation of horton's first law. *Geophysical Bulletin of Hokkaido University* 15 (1966), 1–19. 18

[Wei01] WEISS A.: Topographic position and landforms analysis (poster presentation). In *ESRI International User Conference, July 2001, San Diego, CA* (2001). 4, 9

[Wer88] WERNER C.: Formal analysis of ridge and channel patterns in maturely eroded terrain. *Annals of the Association of American Geographers* 78, 2 (1988), 253–270. 17

[Whe95] WHEATLEY D.: Cumulative viewshed analysis: a gis-based method for investigating intervisibility, and its archaeological application. *Archaeology and Geographic Information Systems: A European perspective* (1995), 171–185. 6

[WL06] WANG L., LIU H.: An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling. *International Journal of Geographical Information Science* 20, 2 (2006), 193–213. 13

[WOB*07] WILSON M. F., O'CONNELL B., BROWN C., GUINAN J. C., GREHAN A. J.: Multiscale terrain analysis of multibeam bathymetry data for habitat mapping on the continental slope. *Marine Geodesy* 30, 1-2 (2007), 3–35. 3, 4

[Woo96] WOOD J.: *The geomorphological characterisation of digital elevation models (PhD thesis)*. University of Leicester (United Kingdom), 1996. 4, 9, 10, 11, 29

[WS78] WISCHMEIER W. H., SMITH D. D.: Predicting rainfall erosion losses: a guide to conservation planning. No. 537. Department of Agriculture, Science and Education Administration, 1978. 15

[WT99] WHIPPLE K. X., TUCKER G. E.: Dynamics of the stream-power river incision model: Implications for height limits of mountain ranges, landscape response timescales, and research needs. *Journal of Geophysical Research: Solid Earth* 104, B8 (1999), 17661–17674. 14

[WYL18] WU J., YANG Q., LI Y.: Partitioning of terrain features based on roughness. *Remote Sensing* 10, 12 (2018), 1985. 5

[Xu22] XU K.: Datumless topography: A universally consistent way to quantify relief, 2022. 16

[YSP02] YOKOYAMA R., SHIRASAWA M., PIKE R. J.: Visualizing topography by openness: a new application of image processing to digital elevation models. *Photogrammetric engineering and remote sensing* 68, 3 (2002), 257–266. 7, 10

[ZOK11] ZAKŠEK K., OŠTIR K., KOKALJ V.: Sky-view factor as a relief visualization technique. *Remote Sensing* 3, 2 (2011), 398–415. 7

[ZSF16] ZHOU G., SUN Z., FU S.: An efficient variant of the priority-flood algorithm for filling depressions in raster digital elevation models. *Computers & Geosciences* 90 (2016), 87–96. 13

[ZSTR07] ZHOU H., SUN J., TURK G., REHG J. M.: Terrain synthesis from digital elevation models. *IEEE Transactions on Visualization and Computer Graphics* 13, 4 (2007), 834–848. 1, 11, 21

## Appendix A: Fractional Laplacian

Let $s \in [0, 1]$, then $\forall u \in S(\mathbb{R}^n)$:

$$(-\Delta)^s u(\mathbf{x}) = C(n, s) \lim_{\epsilon \to 0^+} \int_{\mathbb{R}^n - B(\mathbf{x}, \epsilon)} \frac{u(\mathbf{x}) - u(\mathbf{y})}{\|\mathbf{x} - \mathbf{y}\|^{n+2s}} d\mathbf{y}$$

The normalizing constant $C(n, s)$ is chosen so that different definitions are equivalent. Let $\Gamma$ the Gamma function defined on $\mathbb{R}^+$:

$$\Gamma(z) = \int_{0}^{+\infty} t^{z-1} e^{t} dt$$

We have:

$$C(n, s) = \frac{s4^s \Gamma(s + n/2)}{\pi^{n/2} \Gamma(1 - s)}$$

Moreover, $C(n, s) \equiv s(1-s)$ as $s \to 0^+$ and $s \to 1^-$, and the following identities hold:

$$\lim_{s \to 0^+} (-\Delta)^s u = u \qquad \lim_{s \to 1^-} (-\Delta)^s u = -\Delta u$$

## Appendix B: Quadric fitting

For a symmetric $n \times n$ window, the coefficients are the solution to the following linear system:

$$M_n = \begin{pmatrix} a_{20} \\ a_{02} \\ a_{11} \\ a_{10} \\ a_{01} \\ a_{00} \end{pmatrix} = \begin{pmatrix} \sum z_i x_i^2 \\ \sum z_i y_i^2 \\ \sum z_i x_i y_i \\ \sum z_i x_i \\ \sum z_i y_i \\ \sum z_i \end{pmatrix}$$

$$M_n = \begin{pmatrix} \sum x_i^4 & \sum x_i^2 y_i^2 & 0 & 0 & 0 & \sum x_i^2 \\ \sum x_i^2 y_i^2 & \sum x_i^4 & 0 & 0 & 0 & \sum y_i^2 \\ 0 & 0 & \sum x_i^2 y_i^2 & 0 & 0 & 0 \\ 0 & 0 & 0 & \sum x_i^2 & 0 & 0 \\ 0 & 0 & 0 & 0 & \sum y_i^2 & 0 \\ \sum x_i^2 & \sum y_i^2 & 0 & 0 & 0 & n^2 \end{pmatrix}$$

where $x_i, y_i$ are the local coordinates in the window and $z_i = h(\mathbf{p} + (x_i, y_i))$. Note that $M_n$ only depends on the window offsets, so it can be precomputed for each desired window size $n$. Moreover, we can force the fitted quadric to interpolate at $h(\mathbf{p})$ by setting $a_{00} = 0$ and expressing all heights relative to $h(\mathbf{p})$. We refer the reader to chapter 4.4 in [Woo96] for more details and derivations.