A simple height field. More...

#include <heightfield.h>

Inheritance diagram for HeightField:

Public Member Functions
	HeightField ()
	Empty.

	HeightField (const ScalarField2 &)
	Create a heightfield from a scalar field.

	HeightField (const Box2 &, const QImage &, const double &=0.0, const double &=256.0 *256.0 - 1.0, bool=true)
	Create a heightfield from an image.

	HeightField (const Box2 &, int, int, const double &=0.0)
	Create a flat heightfield.

	HeightField (const Box2 &, int, int, const QVector< double > &)
	Create a heightfield.

	~HeightField ()
	Empty.

double	GetHeight (const Vector2 &) const
	Compute the height at a given position.

Vector	Vertex (const Vector2 &, bool=true) const
	Compute the vertex position on the terrain.

Vector	Vertex (int, int) const
	Compute the vertex corresponding to a given sample.

QVector< Vector >	ArrayVertexes (const QVector< QPoint > &) const
	Compute the coordinates of a set of points on the grid.

Vector	Normal (const Vector2 &, bool=true) const
	Compute the normal for a given position on the terrain.

Vector	Normal (int, int) const
	Compute the normal at a given sample.

double	Slope (int, int) const
	Compute the slope at a given integer point on the terrain.

double	Slope (const QPoint &, const QPoint &) const
	Compute the slope between two points.

double	AverageSlope (int, int) const
	Compute the average slope at a given integer point on the terrain.

double	Aspect (int, int) const
	Compute the aspect (face orientation) at a given integer point on the terrain.

double	K () const
	Compute the Lipschitz constant of the signed distance field defined as f(x,y,z)=z-h(x,y).

void	Subdivide (const double &, Random &=Random::R239)
	Refine the terrain.

ScalarField2	StreamArea (const double &=1.0) const
	Compute the stream area field of the terrain.

ScalarField2	StreamAreaSteepest () const
	Compute the stream area field of the terrain.

ScalarField2	StreamAreaLimited (const double &) const
	Compute the stream area field of the terrain using a limited flow propagation.

ScalarField2	StreamLength (bool=false) const
	Compute the stream length.

ScalarField2	StreamSourceElev (bool=false) const
	Return the maximum source elevation of the flows converging to each point in the height field.

VectorField2	FlowField (ScalarField2 &) const
	Compute the flow field of the terrain.

ScalarField2	AverageSlope () const
	Compute the average slope field of the terrain.

ScalarField2	Slope (bool=false) const
	Compute the maximum slope field of the terrain.

ScalarField2	Aspect () const
	Compute the aspect field (face orientation) of the terrain.

ScalarField2	WetnessIndex (const double &=1.0e-3, bool=false) const
	Compute the wetness index field of the terrain.

ScalarField2	StreamPower (bool=false) const
	Compute the stream power field of the terrain.

ScalarField2	StreamPower (double, double) const
	Compute the stream power field of the terrain.

ScalarField2	ReliefSteepness (const double &) const
	Compute the omni-directional relief and steepness.

ScalarField2	Jut (const double &) const
	Compute the Jut metric.

ScalarField2	Rut (const double &) const
	Compute the Rut metric.

ScalarField2	Peakedness (const double &) const
	Compute the peakedness of a terrain, defined as the percentage of cells in a disk that have lower or equal elevation than the center position.

ScalarField2	TopographicPositionIndex (const double &r) const
	Compute the TPI, an index that measures how much a cell rises with respect to the mean elevation of a given radius. It is very similar in spirit to Peakedness, but here we obtain a value in m instead of the % of cells below it.

ScalarField2	TopographicPositionIndex_SAT (const double &r) const
	Compute the TPI, an index that measures how much a cell rises with respect to the mean elevation of a given radius. It is very similar in spirit to Peakedness, but here we obtain a value in m instead of the % of cells below it.

ScalarField2	TerrainRuggednessIndex (int=3, bool=false) const
	Compute the TRI, an index that measures the ruggedness of a terrain neighborhood by computing the RMS of the residuals with respect to the center location.

ScalarField2	LocalVariance (int=3) const
	Compute the local variance.

ScalarField2	LocalRange (int=3) const
	Compute the local range.

ScalarField2	AreaRatio (int=3) const
	Compute the surface to planar ratio of areas, a measure of terrain ruggedness.

ScalarField2	SurfaceRoughness (int=3) const
	Compute the surface roughness.

ScalarField2	HillslopeAsymmetry (double, double=0, double=45) const
	Compute the hillslope asymmetry map for a given direction.

ScalarField2	Accessibility (const double &, int=16) const
	Compute the accessibility.

ScalarField2	ClearSky (const double &r, int n) const
	Compute the clear sky.

ScalarField2	Light (const Vector &, bool=false) const
	Compute the direct lighting.

int	FillDepressions ()
	Fills all pits and removes all digital dams from the HeightField.

int	FillDepressions (const double &)
	Modifies the HeightField to guarantee drainage.

int	FillDepressions (ScalarField2 &) const
	Fills all pits and removes all digital dams from the HeightField.

int	FillDepressions (const double &, ScalarField2 &) const
	Modifies the HeightField to guarantee drainage.

void	CompleteBreach ()
	Breach depressions.

bool	Intersect (const Ray &, double &, Vector &) const
	Compute the intersection between a ray and the surface of the terrain.

bool	Intersect (const Ray &, double &, Vector &, const Box &, const double &, const double &=1.0e8, const double &=1.0e-4) const
	Compute the intersection between a ray and the surface of the heightfield.

bool	Intersect (const Ray &, QVector< double > &, QVector< Vector > &, const Box &, const double &) const
	Compute the intersection between a ray and the surface of the heightfield.

bool	IntersectBox (const Ray &, double &, Vector &, Vector &, const Box &) const
	Compute the intersection between a ray and the surface border of the heightfield.

double	Shade (const Vector &, const QVector< Vector > &, const QVector< double > &, const double &) const
	Shade a vector of the terrain.

ScalarField2	Sun (const double &, int=5, int=1, int=0, int=12) const
	Compute the average direct lighting from the sun over a year.

ScalarField2	SelfShadow (const Vector &) const
	Compute the self shadowing map.

Array2I	Flow (bool=false) const
	Compute the flow topology.

Array2I	Geomorphons (const double &=0.02) const
	Compute the geomorphon index map.

Array2I	GeomorphonsTangent (double, double=0.02) const
	Compute the geomorphon index map.

Array2I	GeomorphonsAggregated (const double &=0.02) const
	Compute the aggregated geomorphon index map.

Box	GetBox () const
	Get the bounding box of the heightfield.

double	Signed (const Vector &) const
	Compute the signed distance do the heightfield.

void	Scale (const Vector &)
	Scale the height field.

void	Scale (const double &)
	Uniform scaling.

void	Unity ()
	Scale the heightfield so that the compact support should fit within unit box.

void	Draw (QGraphicsScene &, const QPen &=QPen(), const QBrush &=QBrush()) const
	Draw a heightfield.

QVector< QPoint >	Up (int, int) const
	Compute the cells that are uphill of input point.

QVector< QPoint >	Down (int, int) const
	Compute the cells that are downhill of input point.

int	CheckFlowSlope (const QPoint &, FlowStruct &, const double &=1.0) const
	Compute the flow directions at a given point, using an L² metric.

int	CheckFlowDirectionsAngle (const QPoint &, const double &, FlowStruct &) const
	Compute the flow directions at a given point.

Mesh	CreateMesh (bool=true, bool=false, const double &=0.0) const
	Create the geometry of the heightfield.

Mesh	CreateMeshSide (const double &=0.0) const
	Create the geometry of the border of the terrain.

void	CreateCubes (Box &, QVector< FrameScaled > &, const double &=-10.0) const
	Create a stack representation of the model.

void	Carve (const CubicCurve &, const double &, const double &)
	Carve the terrain using a cubic curve, whose elevation define the prescribed elevation along the curve.

void	Carve (const QuadricCurve &, const double &, const double &)
	Carve the terrain using a quadric curve, whose elevation define the prescribed elevation along the curve.

void	Carve (const CubicCurveSet &, const double &, const double &)
	Carve the terrain using a piecewise cubic curve, whose elevation define the prescribed elevation along the curve.

void	Carve (const QuadricCurveSet &, const double &, const double &)
	Carve the terrain using a piecewise quadric curve, whose elevation define the prescribed elevation along the curve.

void	Carve (const Segment &, const double &, const double &)
	Carve the terrain using a segment, whose elevation define the prescribed elevation along the curve.

void	Carve (const SmoothDisc2 &, const double &)
	Flatten the scalar field.

ScalarField2	GradientDistance (const HeightField &) const
	Compute the gradient distance between two heightfields.

ScalarField2	Sea (const double &=0.0) const
	Compute the water elevation in every cell according to a uniform sea level.

QVector< double >	Cross (const Vector2 &, const Vector2 &, int) const
	Compute the elevation along a segment.

void	SmoothBreachMultiScale (const QVector< int > &)
	Smooth breaching with a given input series.

void	SmoothBreachGeometric (double, double)
	Smooth breaching with a geometric series.

void	SmoothBreachLinear (double)
	Smooth breaching process.

void	LowErosionMultiScale (const QVector< double > &, const QVector< double > &)
	Preparametered multi-scale medium scale erosion. Similar capacity limited erosion/deposition to Mei [2007], but with drainage-like flow routing instead of shallow water.

void	LowErosionGeometric (double, double, double, double)
	Preparametered multi-scale medium scale erosion. Similar capacity limited erosion/deposition to Mei [2007], but with drainage-like flow routing instead of shallow water.

QString	EncodeSize () const
	Code the size and range parameters of the heightfield into a string.

Static Public Member Functions
static Box	Size (const QString &)
	Compute the box of the heightfield from a string.

Static Public Attributes
static const double	flat = 1.0e-8
	Small negative epsilon value used in breaching and flow algorithms.

Detailed Description

A simple height field.

Constructor & Destructor Documentation

◆ HeightField() [1/4]

HeightField::HeightField ( const ScalarField2 & s )

Create a heightfield from a scalar field.

This constructor provides implicit conversion.

Parameters

s	Scalar field.

◆ HeightField() [2/4]

HeightField::HeightField	(	const Box2 &	box,
		const QImage &	image,
		const double &	a = 0.0,
		const double &	b = 256.0 * 256.0 - 1.0,
		bool	grayscale = true )

explicit

Create a heightfield from an image.

Parameters

box	Rectangle domain of the terrain.
image	Elevation image.
a,b	Minimum and maximum elevation range.
grayscale	Boolean set to false if the image is provided in color.

◆ HeightField() [3/4]

HeightField::HeightField	(	const Box2 &	box,
		int	nx,
		int	ny,
		const double &	v = 0.0 )

explicit

Create a flat heightfield.

Parameters

box	Rectangle domain of the terrain.
nx,ny	Samples.
v	Constant elevation.

◆ HeightField() [4/4]

HeightField::HeightField	(	const Box2 &	box,
		int	nx,
		int	ny,
		const QVector< double > &	v )

explicit

Create a heightfield.

Parameters

box	Rectangle domain of the terrain.
nx,ny	Samples.
v	Set of elevation values.

Member Function Documentation

◆ Accessibility()

ScalarField2 HeightField::Accessibility	(	const double &	r,
		int	n = 16 ) const

Compute the accessibility.

Parameters

r	Radius.
n	Number of rays.

◆ AreaRatio()

ScalarField2 HeightField::AreaRatio ( int w = 3 ) const

Compute the surface to planar ratio of areas, a measure of terrain ruggedness.

The function computes the ratio of surface triangle areas (formed as 2 triangles every 2x2 cells) with respect to the planar area of those cells. Larger values indicate steeper terrain, minimum value is 1 for flat regions.

Parameters

w	Window size

Author: Oscar Argudo

◆ ArrayVertexes()

QVector< Vector > HeightField::ArrayVertexes ( const QVector< QPoint > & p ) const

Compute the coordinates of a set of points on the grid.

See also: Array2::ArrayVertexes

Parameters

p	Set of point.

◆ Aspect()

double HeightField::Aspect	(	int	i,
		int	j ) const

Compute the aspect (face orientation) at a given integer point on the terrain.

Parameters

i,j	Integer coordinates

◆ AverageSlope()

double HeightField::AverageSlope	(	int	i,
		int	j ) const

Compute the average slope at a given integer point on the terrain.

Simply average the slope in 8 directions.

Parameters

i,j	Integer coordinates

◆ Carve() [1/6]

void HeightField::Carve	(	const CubicCurve &	c,
		const double &	r,
		const double &	b )

Carve the terrain using a cubic curve, whose elevation define the prescribed elevation along the curve.

Parameters

c	Curve.
r	Radius.
b	Blending radius.

◆ Carve() [2/6]

void HeightField::Carve	(	const CubicCurveSet &	c,
		const double &	r,
		const double &	b )

Carve the terrain using a piecewise cubic curve, whose elevation define the prescribed elevation along the curve.

Parameters

c	Curve.
r	Radius.
b	Blending radius.

◆ Carve() [3/6]

void HeightField::Carve	(	const QuadricCurve &	c,
		const double &	r,
		const double &	b )

Carve the terrain using a quadric curve, whose elevation define the prescribed elevation along the curve.

Parameters

c	Curve.
r	Radius.
b	Blending radius.

◆ Carve() [4/6]

void HeightField::Carve	(	const QuadricCurveSet &	c,
		const double &	r,
		const double &	b )

Carve the terrain using a piecewise quadric curve, whose elevation define the prescribed elevation along the curve.

Parameters

c	Curve.
r	Radius.
b	Blending radius.

◆ Carve() [5/6]

void HeightField::Carve	(	const Segment &	c,
		const double &	r,
		const double &	b )

Carve the terrain using a segment, whose elevation define the prescribed elevation along the curve.

Parameters

c	Segment.
r	Radius.
b	Blending radius.

◆ Carve() [6/6]

void HeightField::Carve	(	const SmoothDisc2 &	disc,
		const double &	e )

Flatten the scalar field.

Parameters

disc	Smooth disc reion.
e	Target elevation.

◆ CheckFlowDirectionsAngle()

int HeightField::CheckFlowDirectionsAngle	(	const QPoint &	p,
		const double &	s,
		FlowStruct &	flow ) const

Compute the flow directions at a given point.

Parameters

p	Point.
s	Slope.
flow	Flow information.

◆ CheckFlowSlope()

int HeightField::CheckFlowSlope	(	const QPoint &	p,
		FlowStruct &	flow,
		const double &	power = 1.0 ) const

Compute the flow directions at a given point, using an L² metric.

Parameters

p	Point.
flow	Flow information.
power	Power, set to 1.0 for diffusive flow routing, converges to steepest slope as power increases to infinity.

◆ ClearSky()

ScalarField2 HeightField::ClearSky	(	const double &	r,
		int	n ) const

Compute the clear sky.

Parameters

r	Radius.
n	Number of rays.

◆ CompleteBreach()

void HeightField::CompleteBreach ( )

Breach depressions.

Depression breaching drills a path from a depression's pit cell (its lowest point) along the least-cost (Priority-Flood) path to the nearest cell outside the depression to have the same or lower elevation.

See https://rbarnes.org/

See https://github.com/r-barnes/richdem

See https://github.com/r-barnes/richdem/blob/master/include/richdem/depressions/Lindsay2016.hpp

Author: John Lindsay, implementation by Richard Barnes (rbarn.nosp@m.es@u.nosp@m.mn.ed.nosp@m.u).

◆ CreateCubes()

void HeightField::CreateCubes	(	Box &	cube,
		QVector< FrameScaled > &	frames,
		const double &	z = -10.0 ) const

Create a stack representation of the model.

Stacks are defined as scaled boxes.

Parameters

cube	Reference cube.
frames	Set of instances.
z	Base elevation.

◆ CreateMesh()

Mesh HeightField::CreateMesh	(	bool	surface = true,
		bool	side = false,
		const double &	z = 0.0 ) const

Create the geometry of the heightfield.

Parameters

surface	Generate the surface elevation of the heightfield.
side	Generate the side of the heightfield.
z	The height of the base of the block.

◆ CreateMeshSide()

Mesh HeightField::CreateMeshSide ( const double & z = 0.0 ) const

Create the geometry of the border of the terrain.

Parameters

z	The height of the base of the block.

◆ Cross()

QVector< double > HeightField::Cross	(	const Vector2 &	a,
		const Vector2 &	b,
		int	n ) const

Compute the elevation along a segment.

Parameters

a,b	Segment.
n	Number of samples.

◆ Down()

QVector< QPoint > HeightField::Down	(	int	x,
		int	y ) const

Compute the cells that are downhill of input point.

Parameters

x,y	Integer coordinates of the point.

◆ Draw()

void HeightField::Draw	(	QGraphicsScene &	scene,
		const QPen &	pen = QPen(),
		const QBrush &	brush = QBrush() ) const

Draw a heightfield.

Parameters

scene	Graphics scene.
pen	The pen.
brush	The brush.

◆ EncodeSize()

QString HeightField::EncodeSize ( ) const

Code the size and range parameters of the heightfield into a string.

output / encode a string : x=40800 y=40800 z=+127 +3687 in meters or c=25645-z= in centimeters for cell size and grey level

◆ FillDepressions() [1/4]

int HeightField::FillDepressions ( )

Fills all pits and removes all digital dams from the HeightField.

See also: FillDepressions(ScalarField2&) const

Author: Mathieu Gaillard

Returns: The number of cells where a significant alteration of the DEM has occurred.

◆ FillDepressions() [2/4]

int HeightField::FillDepressions ( const double & eps )

Modifies the HeightField to guarantee drainage.

See also: FillDepressions(const double&, ScalarField2&) const

Author: Mathieu Gaillard

Parameters

eps	The minimal elevation difference between two cells.

Returns: The number of cells where a significant alteration of the DEM has occurred.

◆ FillDepressions() [3/4]

int HeightField::FillDepressions	(	const double &	eps,
		ScalarField2 &	depression_field ) const

Modifies the HeightField to guarantee drainage.

Algorithm: Priority-flood: An optimal depression-filling and watershed-labeling algorithm for digital elevation models. Barnes, R., Lehman, C., Mulla, D., 2014.

This version of Priority-Flood starts on the edges of the DEM and then works its way inwards using a priority queue to determine the lowest cell which has a path to the edge. The neighbours of this cell are added to the priority queue if they are higher. If they are lower, then their elevation is increased by a small amount to ensure that they have a drainage path and they are added to a "pit" queue which is used to flood pits. Cells which are higher than a pit being filled are added to the priority queue. In this way, pits are filled without incurring the expense of the priority queue.

The values to fill the depressions are added to the depression field.

The depression field should have the exact same size as the HeightField. If the given depression field does not have the exact same size as the height field, returns the total number of vertices in the heightfield.

To be sure that the height field is never flat, elevations are always increased by std::nextafter(elevation) + eps. Thus, even if eps == 0, the terrain is never flat.

Warning, because the method is const it is not recommended to use it like that

// Not recommended, although it works.

hf.FillDepressions(0.0, hf);

Author: Mathieu Gaillard

Parameters

eps	The minimal elevation difference between two cells.
depression_field	A scalar field in which the values to fill depressions are added

Returns: The number of cells where a significant alteration of the DEM has occurred.

◆ FillDepressions() [4/4]

int HeightField::FillDepressions ( ScalarField2 & d ) const

Fills all pits and removes all digital dams from the HeightField.

Algorithm: Priority-flood: An optimal depression-filling and watershed-labeling algorithm for digital elevation models. Barnes, R., Lehman, C., Mulla, D., 2014.

Priority-Flood starts on the edges of the DEM and then works its way inwards using a priority queue to determine the lowest cell which has a path to the edge. The neighbours of this cell are added to the priority queue if they are higher. If they are lower, they are raised to the elevation of the cell adding them, thereby filling in pits. The neighbors are then added to a "pit" queue which is used to flood pits. Cells which are higher than a pit being filled are added to the priority queue. In this way, pits are filled without incurring the expense of the priority queue.

The values to fill the depressions are added to the depression field.

The depression field should have the exact same size as the HeightField. If the given depression field does not have the exact same size as the height field, returns the total number of vertices in the heightfield.

Warning, because the method is const it is not recommended to use it like that

// Not recommended, although it works.

hf.FillDepressions(hf);

Author: Mathieu Gaillard

Parameters

d	A scalar field in which the values to fill depressions are added

Returns: The number of cells where a significant alteration of the DEM has occurred.

◆ Flow()

Array2I HeightField::Flow ( bool steep = false ) const

Compute the flow topology.

Parameters

steep Flag, set to true if we rely on the steepest slope, false for multiple convergent and divergent flows.

◆ FlowField()

VectorField2 HeightField::FlowField ( ScalarField2 & stream ) const

Compute the flow field of the terrain.

The direction if the average downhill flow.

This function differs slightly from HeightField::StreamArea(). Pits have a null flow vector, although the stream area is a local maximum.

The stream area is computed as a fraction of cells, independtly of the surface area of a cell.

Parameters

stream Returned stream area;

◆ Geomorphons()

Array2I HeightField::Geomorphons ( const double & t = 0.02 ) const

Compute the geomorphon index map.

Geomorphon classify cells into ten categories: 0: flat, 1: peak, 2: ridge, 3: shoulder, 4:hollow=convex, 5: slope, 6: spur=concave, 7: footslope, 8: valley, and 9: pit.

Parameters

t	Slope threshold.

◆ GeomorphonsAggregated()

Array2I HeightField::GeomorphonsAggregated ( const double & t = 0.02 ) const

Compute the aggregated geomorphon index map.

Aggregation classifies into only five categories: flat (flat), ridges (ridge, peak, spur, shoulder), slope (slope), valley (valley, pit, hollow) and foot slope (foot slope).

Parameters

t	Slope threshold.

◆ GeomorphonsTangent()

Array2I HeightField::GeomorphonsTangent	(	double	maxDist,
		double	flatTangent = 0.02 ) const

Compute the geomorphon index map.

Geomorphon classify cells into ten categories: flat, peak, ridge, shoulder, hollow, slope, spur, footslope, valley, and pit.

Parameters

maxDist	maximum ray traversal distance.
flatTangent	threshold for considering a ray angle as flat.

◆ GetBox()

Box HeightField::GetBox ( ) const

Get the bounding box of the heightfield.

Although this function has the same name as Array2::GetBox(), it computes the minimum and maximum elevation of the terrain (computationally intensive).

The two-dimensional box can be obtained by using:

HeightField heightfield;

Box2 box = heightfield.Array2::GetBox(); // Get the domain in the plane.

Box2::Box2

Box2()

HeightField::HeightField

HeightField()

Empty.

Definition heightfield.h:105

See also: Array2::GetBox()

◆ GetHeight()

double HeightField::GetHeight ( const Vector2 & p ) const

inline

Compute the height at a given position.

Parameters

p Point.

◆ GradientDistance()

ScalarField2 HeightField::GradientDistance ( const HeightField & h ) const

Compute the gradient distance between two heightfields.

Parameters

h	Second heightfield.

◆ HillslopeAsymmetry()

ScalarField2 HeightField::HillslopeAsymmetry	(	double	km,
		double	direction = 0,
		double	tolerance = 45 ) const

Compute the hillslope asymmetry map for a given direction.

Poulos et al 2012. Hillslope asymmetry maps reveal widespread, multi-scale organization

The authors define hillsope asymmetry as the difference in median slope between slopes of opposite orientation inside an area.

NOTE: I used average slopes instead of median. First, for a reason of efficiency. Second, it also makes more sense in my opinion. In any case, results visually very similar. It is possible to change the behavior in the code, just set MEDIAN_ASYMMETRY = true

Parameters

km	Side length of area in which analysis is performed.
direction	Angle (in degrees, 0 is +X).
tolerance	How much we can deviate from this angle.

Author: Oscar Argudo

◆ Intersect() [1/3]

bool HeightField::Intersect	(	const Ray &	ray,
		double &	t,
		Vector &	q ) const

Compute the intersection between a ray and the surface of the terrain.

The algorithm uses a ray marching approach, therefore this function may require many iterations and the resulting intersection may not be accurate.

Parameters

ray	The ray.
t	Intersection depth.
q	Intersection point.

◆ Intersect() [2/3]

bool HeightField::Intersect	(	const Ray &	ray,
		double &	t,
		Vector &	q,
		const Box &	box,
		const double &	k,
		const double &	length = 1.0e8,
		const double &	epsilon = 1.0e-4 ) const

Compute the intersection between a ray and the surface of the heightfield.

The algorithm uses a ray marching approach, therefore this function may require many iterations and the resulting intersection may not be accurate.

Parameters

ray	The ray.
t	Returned distance along the ray.
box	Box where intersection will be computed.
q	Returned intersection point.
k	Lipschitz constant.
length	Maximum distance along the ray.
epsilon	Minimum stepping distance.

◆ Intersect() [3/3]

bool HeightField::Intersect	(	const Ray &	ray,
		QVector< double > &	tv,
		QVector< Vector > &	q,
		const Box &	box,
		const double &	k ) const

Compute the intersection between a ray and the surface of the heightfield.

The algorithm uses a sphere tracing approach, therefore this function may require many iterations and the resulting intersection may not be accurate.

Parameters

ray	The ray.
tv	Returned distance along the ray.
box	Box where intersection will be computed.
q	Returned intersection point.
k	Lipschitz constant.

◆ IntersectBox()

bool HeightField::IntersectBox	(	const Ray &	ray,
		double &	t,
		Vector &	p,
		Vector &	n,
		const Box &	box ) const

Compute the intersection between a ray and the surface border of the heightfield.

Parameters

ray	The ray.
t	Distance along the ray.
box	Box where intersection will be computed.
p	Returned intersection point.
n	Returned normal.

◆ Jut()

ScalarField2 HeightField::Jut ( const double & r ) const

Compute the Jut metric.

See Xu 2022, "Datumless Topography, A Universally Consistent Way to Quantify Relief".

It is defined as the maximum angle-reduced height of p with respect to any other point q. Jut reaches larger values - more impressiveness - when the terrain rises abruptly in a small area, such as in the vicinity of cliffs. Note that this metric is intended to be used taking into account planetary curvature, here we coded a simplified version assuming the height field domain is a plane.

Parameters

r	Radius for cutting the computation at some point.

Author: Oscar Argudo

◆ Light()

ScalarField2 HeightField::Light	(	const Vector &	u,
		bool	cosine = false ) const

Compute the direct lighting.

Parameters

u	Light direction.
cosine	Boolean, set to true to use scalar product between normal and light, false to use (1+c)/2.

◆ LocalRange()

ScalarField2 HeightField::LocalRange ( int w = 3 ) const

Compute the local range.

Range is max - min inside a window.

Parameters

w	Window size (authors recommend 3).

Author: Oscar Argudo

◆ LocalVariance()

ScalarField2 HeightField::LocalVariance ( int w = 3 ) const

Compute the local variance.

Woodcock and Strahler 1987. The factor of scale in remote sensing. Dragut et al. 2011. Local variance for multi-scale analysis in geomorphometry.

Local variance is defined as the standard deviation of a window centered around the cell.

Compute the local variance of every cell and take the average over all the domain.

Authors plot the mean local variance with respect to image scale, and show that peaks in this graph occur at 1/2 to 3/4 the size of objects in the image, so it could be used as a scale indicator.

Parameters

w	Window size (authors recommend 3).

Author: Oscar Argudo

◆ LowErosionGeometric()

void HeightField::LowErosionGeometric	(	double	max_radius,
		double	radius_attenuation,
		double	max_iter,
		double	iter_attenuation )

Preparametered multi-scale medium scale erosion. Similar capacity limited erosion/deposition to Mei [2007], but with drainage-like flow routing instead of shallow water.

Parameters

max_radius	Largest erosion radius.
radius_attenuation	Attenuation coeff to compute next radii.
max_iter	Number of iteration for the largest radius.
iter_attenuation	Attenuation coeff to compute next number of iteration. Stops when radius is one pixel or below.

◆ LowErosionMultiScale()

void HeightField::LowErosionMultiScale	(	const QVector< double > &	radii,
		const QVector< double > &	iters )

Preparametered multi-scale medium scale erosion. Similar capacity limited erosion/deposition to Mei [2007], but with drainage-like flow routing instead of shallow water.

Parameters

radii	List of erosion radius for each scale.
iters	List of number of iterations to apply at each scale.

◆ Normal() [1/2]

Vector HeightField::Normal	(	const Vector2 &	p,
		bool	triangular = true ) const

Compute the normal for a given position on the terrain.

Note that this function may be expensive to compute.

Parameters

p	Point.
triangular	Boolean, use triangle normals if set to true, bilinear interpolation of normals at vertices otherwize.

◆ Normal() [2/2]

Vector HeightField::Normal	(	int	i,
		int	j ) const

Compute the normal at a given sample.

This function uses the weighted sum (area) of the normals of the triangles sharing the point on the grid. The returned vector is normalized.

Parameters

i,j	Integer coordinates of the sample.

◆ Peakedness()

ScalarField2 HeightField::Peakedness ( const double & r ) const

Compute the peakedness of a terrain, defined as the percentage of cells in a disk that have lower or equal elevation than the center position.

See Llobera 2001, Building Past Landscape Perception With GIS: Understanding Topographic Prominence. Although the author calls it prominence, it is not the same definition as used in orometry for peak prominence.

See Arge et al. 2013: Algorithms for Computing Prominence on Grid Terrains for faster implementations.

Parameters

r	Radius of the disk (we actually use a squared window).

Author: Oscar Argudo

◆ ReliefSteepness()

ScalarField2 HeightField::ReliefSteepness ( const double & r ) const

Compute the omni-directional relief and steepness.

See Earl and Metzler 2015, Cloud-Capped Towers: Capturing Terrain Characteristics Using Topographic Functionals.

From a theoretical point of view, ORS is the value that results from convergence when we integrate all distances. ORS monotonically increases as we make this cut-off bigger, since the function we integrate is always positive. From my experience (and by definition) after only a few km of distance 10-20, the value is already a good approximation but also, different radius cut-off highlight different interesting properties. If we use a small radius, those points with a very steep edge nearby or even needle-shaped secondary small peaks will have large values of ORS, more than other bigger but gentler peaks. When the radius increases, bigger peaks shadow them with larger ORS values.

Parameters

r	Radius for cutting the computation at some point.

Author: Oscar Argudo

◆ Rut()

ScalarField2 HeightField::Rut ( const double & r ) const

Compute the Rut metric.

See Xu 2022, "Datumless Topography, A Universally Consistent Way to Quantify Relief".

It is defined as the maximum angle-reduced height of any point q with respect to p. Rut measures how sharply or impressively a point dips below its immediate surroundings. Note that this metric is intended to be used taking into account planetary curvature, here we coded a simplified version assuming the height field domain is a plane.

Parameters

r	Radius for cutting the computation at some point.

Author: Oscar Argudo

◆ Scale() [1/2]

void HeightField::Scale ( const double & s )

Uniform scaling.

Parameters

s	Scaling factor.

◆ Scale() [2/2]

void HeightField::Scale ( const Vector & s )

Scale the height field.

The domain is scaled using the x and y components of the scaling vector, whereas the heights of the heightfield are scaled using the z component.

Parameters

s	Scaling factor.

◆ Sea()

ScalarField2 HeightField::Sea ( const double & height = 0.0 ) const

Compute the water elevation in every cell according to a uniform sea level.

Returns: The amount of water in every cell.

◆ SelfShadow()

ScalarField2 HeightField::SelfShadow ( const Vector & light ) const

Compute the self shadowing map.

Parameters

light Light direction.

◆ Shade()

double HeightField::Shade	(	const Vector &	p,
		const QVector< Vector > &	lights,
		const QVector< double > &	intensity,
		const double &	sum ) const

Shade a vector of the terrain.

Parameters

p	Terrain coordinates.
lights	Set of lights.
intensity	Set of light intensities.
sum	Sum.

◆ Signed()

double HeightField::Signed ( const Vector & p ) const

Compute the signed distance do the heightfield.

This function uses ScalarField2::BiCubicValue to compute the elevation.

The volume delimited by the heightfield is embedded an infinite vertical box.

Parameters

p Point.

◆ Size()

Box HeightField::Size ( const QString & s )

static

Compute the box of the heightfield from a string.

Parameters

s String.

◆ Slope() [1/3]

ScalarField2 HeightField::Slope ( bool a = false ) const

Compute the maximum slope field of the terrain.

Simply compute the norm of the gradient.

This function may be used to perform hill-slope erosion simulation.

Parameters

a	Boolean, use average slope if set to true, and slope otherwize.

See also: (), Slope(), AverageSlope(), GradientNorm()

◆ Slope() [2/3]

double HeightField::Slope	(	const QPoint &	pa,
		const QPoint &	pb ) const

Compute the slope between two points.

Parameters

pa,pb Points.

◆ Slope() [3/3]

double HeightField::Slope	(	int	i,
		int	j ) const

Compute the slope at a given integer point on the terrain.

This is a convenience function corresponding to the following piece of code:

double s=Norm(heightfield.Gradient(i, j));

Parameters

i,j	Integer coordinates

◆ SmoothBreachGeometric()

void HeightField::SmoothBreachGeometric	(	double	r,
		double	a )

Smooth breaching with a geometric series.

Parameters

r	Radius, breaching will be performed at the corresponding resolution.
a	Attenuation.

Author: Hugo Schott

◆ SmoothBreachLinear()

void HeightField::SmoothBreachLinear ( double r )

Smooth breaching process.

Parameters

r	Radius, breaching will be performed at the corresponding resolution with linear attenuation.

Author: Hugo Schott

◆ SmoothBreachMultiScale()

void HeightField::SmoothBreachMultiScale ( const QVector< int > & ni )

Smooth breaching with a given input series.

Parameters

ni	Set of blurring iterations.

Author: Hugo Schott

◆ StreamArea()

ScalarField2 HeightField::StreamArea ( const double & power = 1.0 ) const

Compute the stream area field of the terrain.

The stream area is computed as a fraction of cells, independtly of the surface area of a cell.

Parameters

power Power.

◆ StreamAreaLimited()

ScalarField2 HeightField::StreamAreaLimited ( const double & r ) const

Compute the stream area field of the terrain using a limited flow propagation.

Parameters

r	Flow propagation range, this function is the same as HeightField::StreamArea() if range is infinitely large.

See also: HeightField::StreamArea()

◆ StreamAreaSteepest()

ScalarField2 HeightField::StreamAreaSteepest ( ) const

Compute the stream area field of the terrain.

See also: HeightField::StreamArea()

◆ StreamLength()

ScalarField2 HeightField::StreamLength ( bool steepest = false ) const

Compute the stream length.

The stream length is computed, for all cells, as the maximum distance to a spring, which are the peaks.

◆ StreamPower() [1/2]

ScalarField2 HeightField::StreamPower ( bool a = false ) const

Compute the stream power field of the terrain.

The algorithm works as follows: first compute the stream area for every terrain sample, and then compute the stream power index.

Parameters

a	Boolean, use average slope if set to true, and slope otherwize.

◆ StreamPower() [2/2]

ScalarField2 HeightField::StreamPower	(	double	m,
		double	n ) const

Compute the stream power field of the terrain.

The algorithm works as follows: first compute the stream area for every terrain sample, and then compute the stream power index.

Parameters

m,n	Exponent for the stream area and slope, n should be twice as m.

◆ StreamSourceElev()

ScalarField2 HeightField::StreamSourceElev ( bool steepest = false ) const

Return the maximum source elevation of the flows converging to each point in the height field.

Parameters

steepest If true, uses D8 routing algorithm. Otherwise, MFD

◆ Subdivide()

void HeightField::Subdivide	(	const double &	e,
		Random &	r = Random::R239 )

Refine the terrain.

This function subdivides the samples and uses a smooth interpolation. A random elevation may be added to the samples.

Parameters

e	Amplitude of the random elevation added to the samples.
r	Random number generator.

◆ Sun()

ScalarField2 HeightField::Sun	(	const double &	latitude,
		int	daystep = 5,
		int	hourstep = 1,
		int	dayoffset = 0,
		int	houroffset = 12 ) const

Compute the average direct lighting from the sun over a year.

Parameters

latitude	Latitude of the terrain, longitude is set to 0.0.
daystep,hourstep	Day and hour steps.
dayoffset,houroffset	Starting day and hour.

◆ SurfaceRoughness()

ScalarField2 HeightField::SurfaceRoughness ( int w = 3 ) const

Compute the surface roughness.

Lindsay et al 2019. Scale-optimized surface roughness for topographic analysis.

The authors associate surface normals dispersion to roughness, a measure of texture complexity. In particular, they measure the spherical standard deviation in a window w × w, i.e. the angular spread of the normal directions in the window. It will be 0 for flat and inclined planes, and increase with surface texture complexity.

They also propose to remove lower frequencies by doing a gaussian low-pass filter on the DEM, to obtain a signature that depends only on current scale (otherwise, it increases monotonically with w).

Parameters

w	Window size.

Author: Oscar Argudo

◆ TerrainRuggednessIndex()

ScalarField2 HeightField::TerrainRuggednessIndex	(	int	w = 3,
		bool	absolute = false ) const

Compute the TRI, an index that measures the ruggedness of a terrain neighborhood by computing the RMS of the residuals with respect to the center location.

See Riley et al. 1999: A terrain ruggedness index that quantifies topographic heterogeneity

Parameters

w	Window size (authors used 3).
absolute	If true, instead of root of squared sums does the average of absolute value diffs.

Author: Oscar Argudo

◆ TopographicPositionIndex()

ScalarField2 HeightField::TopographicPositionIndex ( const double & r ) const

Compute the TPI, an index that measures how much a cell rises with respect to the mean elevation of a given radius. It is very similar in spirit to Peakedness, but here we obtain a value in m instead of the % of cells below it.

This function is widely used (or was) to classify landforms, and it is integrated in many GIS packages.

See Weiss 2001: Topographic position and landforms analysis (poster) And Wilson and Gallant 2000: Primary Topographic Attributes (chapter 3 in "Terrain Analysis: Principles and Applications")

Parameters

r	Radius of the disk (we actually use a squared window).

Author: Oscar Argudo

◆ TopographicPositionIndex_SAT()

ScalarField2 HeightField::TopographicPositionIndex_SAT ( const double & r ) const

Compute the TPI, an index that measures how much a cell rises with respect to the mean elevation of a given radius. It is very similar in spirit to Peakedness, but here we obtain a value in m instead of the % of cells below it.

Optimized version for large w using a Summed Area Table.

Parameters

r	Radius of the disk (we actually use a squared window).

Author: Oscar Argudo

◆ Unity()

void HeightField::Unity ( )

Scale the heightfield so that the compact support should fit within unit box.

Centers the terrain and scales it.

◆ Up()

QVector< QPoint > HeightField::Up	(	int	x,
		int	y ) const

Compute the cells that are uphill of input point.

Parameters

x,y	Integer coordinates of the point.

◆ Vertex() [1/2]

Vector HeightField::Vertex	(	const Vector2 &	p,
		bool	triangular = true ) const

Compute the vertex position on the terrain.

Parameters

p	Point.
triangular	Boolean, use triangular interpolation if set to true, bilinear interpolation otherwize.

◆ Vertex() [2/2]

Vector HeightField::Vertex	(	int	i,
		int	j ) const

Compute the vertex corresponding to a given sample.

Parameters

i,j	Integer coordinates of the sample.

◆ WetnessIndex()

ScalarField2 HeightField::WetnessIndex	(	const double &	c = 1.0e-3,
		bool	a = false ) const

Compute the wetness index field of the terrain.

The algorithm works as follows: first compute the stream area for every terrain sample, and then compute the wetness index.

Note that the wetness index is defined as ln ( A / s ) where s is the slope, however the definition does not work well for small slopes, therefore the implementation defines the wetness index as ln ( A / ( e + s ) ), where e is 0.001.

Parameters

c	Coefficent for slope factor, if set to 0 the function will compute the logarithm of the drainage area.
a	Boolean, use average slope if set to true, and slope otherwize.

See also: Slope(), AverageSlope()

Public Member Functions

Static Public Member Functions

Static Public Attributes

Detailed Description

Constructor & Destructor Documentation

◆ HeightField() [1/4]

◆ HeightField() [2/4]

◆ HeightField() [3/4]

◆ HeightField() [4/4]

Member Function Documentation

◆ Accessibility()

◆ AreaRatio()

◆ ArrayVertexes()

◆ Aspect()

◆ AverageSlope()

◆ Carve() [1/6]

◆ Carve() [2/6]

◆ Carve() [3/6]

◆ Carve() [4/6]

◆ Carve() [5/6]

◆ Carve() [6/6]

◆ CheckFlowDirectionsAngle()

◆ CheckFlowSlope()

◆ ClearSky()

◆ CompleteBreach()

◆ CreateCubes()

◆ CreateMesh()

◆ CreateMeshSide()

◆ Cross()

◆ Down()

◆ Draw()

◆ EncodeSize()

◆ FillDepressions() [1/4]

◆ FillDepressions() [2/4]

◆ FillDepressions() [3/4]

◆ FillDepressions() [4/4]

◆ Flow()

◆ FlowField()

◆ Geomorphons()

◆ GeomorphonsAggregated()

◆ GeomorphonsTangent()

◆ GetBox()

◆ GetHeight()

◆ GradientDistance()

◆ HillslopeAsymmetry()

◆ Intersect() [1/3]

◆ Intersect() [2/3]

◆ Intersect() [3/3]

◆ IntersectBox()

◆ Jut()

◆ Light()

◆ LocalRange()

◆ LocalVariance()

◆ LowErosionGeometric()

◆ LowErosionMultiScale()

◆ Normal() [1/2]

◆ Normal() [2/2]

◆ Peakedness()

◆ ReliefSteepness()

◆ Rut()

◆ Scale() [1/2]

◆ Scale() [2/2]

◆ Sea()

◆ SelfShadow()

◆ Shade()

◆ Signed()

◆ Size()

◆ Slope() [1/3]

◆ Slope() [2/3]

◆ Slope() [3/3]

◆ SmoothBreachGeometric()

◆ SmoothBreachLinear()

◆ SmoothBreachMultiScale()

◆ StreamArea()

◆ StreamAreaLimited()

◆ StreamAreaSteepest()

◆ StreamLength()

◆ StreamPower() [1/2]

◆ StreamPower() [2/2]

◆ StreamSourceElev()