DEM visualisation techniques: Multi-scale integral invariants

Earlier this year, I was impressed by Hubert Mara‘s presentation at the CAA workshop in Berlin. He had used a method called multi-scale integral invariants (MSII) to extract the incised charcaters from cuneiform tablets and inscription from old tombstones. This was surely be something that could be useful for visualisaing LIDAR-based DEMs: back in the office, I implemented it as an additinal algorithm in LiVT.

Now, how does it work? To begin with, it is a multi-scale approach (hence the name). The algorithm places multiple spheres of different diameters on each pixel in the DEM and computes how much of the volume of the spheres is above/below the DEM surface. As a result, you get a number of values (volume fractions above surface) for each DEM pixel. These sets of n values are interpreted as n-dimensional vectors.

By computing the distance of these n-dimensional vectors from a reference vector, the data can be reduced to a raster map containing a single value for each pixel. Low values (low vector distance) indicate high similarity with the reference vector and vice versa. Using an appropriate greyscale histogram stretch, this raster map can be displayed as an image. The reference vector can, for example, be determined by extracting the vector values for a specific relief feature or a point within a cuneiform character or simply by chosing the origin of the n-dimensional coordinate system (i.e. zero).

MSII visualisation (distance from origin, 8 scales starting with radius=2 and scale-to-scale factor 1.414) of the same area as in the post about local dominance  LIDAR data (c) LGL/LAD.

MSII visualisation (distance from origin, 8 scales starting with radius=2 and scale-to-scale factor 1.414; greyscale histogram stretch 1.35…1.55) of the same area as in the post about local dominance. LIDAR data (c) LGL/LAD.

References

Mara, H., Krömker, S., Jakob, S., Breuckmann, B., 2010. GigaMesh and Gilgamesh – 3D Multiscale Integral Invariant Cuneiform Character Extraction, in: A. Artusi, M. Joly-Parvex, G. Lucet, A. Ribes u. D. Pitzalis (Hg.), The 11th International Symposium on Virtual Reality, Archaeology and Cultural Heritage VAST (Paris, France, 2010), 131–138.

Advertisements

DEM visualisation techniques: Local dominance

Local dominance visualisation of a DEM is based on computing, for every pixel of the DEM, how dominant an observer standing on that point would be for a local surrounding area. Dominance as used here is the average steepness of the angle at which the observer looks down at the surrounding land surface. It is higher for points on local elevations and lower for points in local depressions.

Local dominance is computed for pixels within a specified maximum radius and a specified observer height above the surface. To reduce the noisy appearance of the resulting image due to small-scale surface roughness, a minimum radius can be specified. Pixel brightness is derived from the local dominance values by applying an appropriate greyscale histogram stretch.

Local dominance visualisation is well suited for very subtle positive relief features such as former field boundaries or strongly eroded burial mounds, but also delivers very good results for topographic depressions such as dolines, mining traces or hollow ways.

Shaded relief image (illumination azimuth 45°, elevation 30°, no vertical exaggeration) of a low-relief area. Grey and black lines are present-day roads and field boundaries, respectively. Former field boundaries are faintly visible depending on illumination direction.

Shaded relief image (illumination azimuth 45°, elevation 30°, no vertical exaggeration) of a low-relief area. Grey and black lines are present-day roads and field boundaries, respectively. Former field boundaries are faintly visible depending on illimation direction. LIDAR data (c) LGL/LAD.

Local dominance visualisation of the same area. Former field boundaries are clearly visible.

Local dominance visulaisation of the same area. Former field boundaries and a former road are clearly visible. LIDAR data (c) LGL/LAD.

DEM visualisation techniques: Cumulative visibility

A viewshed is the area visible from a given vantage point. The viewshed area depends on the topographic position of the vantage point and the surrounding topography, but also on the height of the observer standing on the vantage point, the height of the objects that should be visible to the observer (e.g. the ground surface or a second person) and the radius under consideration.

Cumulative visibility, on the other hand, specifies the size of the area from which a point in the DEM (or an object on that point) is visible to observers of a certain height. Essentially, it’s the same principle of intervisibility, just don’t get observer and object height mixed up.

DEM visualisation by cumulative visibility is based on computing, for each pixel of the DEM, the size of the area within a specified radius from which an object is visible. Because the surrounding topography plays a dominant role for intervisibility, the resulting raster map can also be a suitable technique to visualise that topography. Besides this, such a visualisation can be used as a tool for analysing for example the locations of archaeological sites.

The resulting raster map contains percentage values (0…100) of the size of the cumulative visibility area relative to the entire area within the specified radius.

Visualisation of cumulative visibility (radius 100 m, greyscale histogram stretch to 10...75%) of the same area as in the previous post on accessibility visualisation.  LIDAR data (c) LGL/LAD.

Visualisation of cumulative visibility (radius 100 m, greyscale histogram stretch to 10…75%) of the same area as in the previous post on accessibility visualisation. LIDAR data (c) LGL/LAD.

DEM visualisation techniques: Accessibility

DEM data can be visualised by computing surface accessibility. This means that an algorithm determines (for every pixel of the DEM) the maximum radius of a sphere that could be placed on the surface at this position without being impeded by the heights of surrounding pixels.

To reduce computation time, the algorithm only takes into account surrounding pixels within a pre-defined radius. Computation time can be further reduced by taking into account only surrounding pixels along a small number of radial lines rather than all pixels.

The range of values in the resulting accessibility raster map corresponds to the range of sphere radii. Greycale or colour mapping is used to display the results as an image.

Shaded relief image (ilumination azimuth 315°, elevation 45°, no vertical exaggeration) of an area on the edge of the Upper Rhine Valley. LIDAR data (c) LGL/LAD.

Shaded relief image (ilumination azimuth 315°, elevation 45°, no vertical exaggeration) of an area on the edge of the Upper Rhine Valley. LIDAR data (c) LGL/LAD.

Accessibility visualisation (greyscale histogram stretch for accessibility radii 1...15 m) of the same area. LIDAR data (c) LGL/LAD.

Accessibility visualisation (greyscale histogram stretch for accessibility radii 1…15 m) of the same area. LIDAR data (c) LGL/LAD.

Accessibility is particularly suitable for visualising negative relief features (e.g. pits, hollow ways) and features on slopes (e.g. agricultural terraces). Subtle relief features on more or less horizontal surfaces show up only poorly or not at all.

As an additional modification of the algorithm, the implementation in the Lidar Visualisation Toolbox LiVT also allows the computation of surface accessibility for directional cones.

References

Miller, G., 1994. Efficient algorithm for local and global accessibility shading. Computer Graphics Proceedings, Annual Conference Series SIGGRAPH, 319–325.

DEM visualisation techniques: Local relief model

Trend removal is a very useful tool to extract topographic detail from DEMs. It can help us to distinguish between positive and negative deviations or “anomalies” from the general (smoothed) forms of the earth surface. But what about measuring the height or depth of these deviations? For trend removal, we simply subtract a smoothed version of the DEM from the original DEM. Any local topographic “anomalies” will smoothed rather than being removed completely. Because of this, more prominent relief features will to a certain extent still be present in smoothed version of the DEM. Ad a result of this, the relative height of such local relief features will be underestimated.

To be able to display (and/or measure) the relative height/depth of local relief “anomalies” it would therefore be good to cut out these features from the DEM and then subtract this “purged” DEM from the original DEM. Trend removal is used as a first step in this process: it results in a difference map containing positive and negative values for positive and negative local relief features, respectively.

In a second step, we can delineate the boundaries between positive and negative local relief features by finding the positions where trend removal values are zero. If we extract the DEM elevation values at these positions and fill the interjacent areas using interpolation, we have a purged DEM in which the local relief features are actually cut out.

Shaded relief image (illumination azimuth 315°, elevation 45°, no vertical exaggeration) overlaid with the zero contours derived from the trend removal map.

Shaded relief image (illumination azimuth 315°, elevation 45°, no vertical exaggeration) overlaid with the zero contours derived from the trend removal map. LIDAR data (c) LGL/LAD.

Shaded relief image (illumination azimuth 315°, elevation 45°, no vertical exaggeration) of the purged DEM. Note that while almost all relief detail has been cut out, part of the largest burial mound is still present in the purged DEM because its diameter  strongly exceeds the ow-pass filter diameter used for trend removal.

Shaded relief image (illumination azimuth 315°, elevation 45°, no vertical exaggeration) of the purged DEM. Note that while almost all relief detail has been cut out, part of the largest burial mound is still present in the purged DEM because its diameter strongly exceeds the ow-pass filter diameter used for trend removal. LIDAR data (c) LGL/LAD.

By subtracting this purged DEM from the original DEM we get what I call a Local Relief Model (LRM), which is a digital model of local relief features. This can be displayed as an image using greyscale or colour mapping. I usually use blue colours for negative and red to yellow colours for positive relief features to make it easier to intuitively read the image.

Colour-coded LRM. Note the strong circular artefact around the largest burial mound. The profile across a smaller burial mound conveys realistic relative height values.

Colour-coded LRM. Note the strong circular artefact around the largest burial mound. The profile across a smaller burial mound conveys realistic relative height values. LIDAR data (c) LGL/LAD.

The LRM does represent the height and depth of local relief features more realistically than the simple trend removal. However, because it is based on the same general approach it also suffers from the same computation artefacts as simple trend removal: there can be apparent ditches or banks which are not actually present.

When working with LRM, I usually display a colour-coded LRM as an overlay over (i.e., multiply by) a shaded relief image. That way, I can visually combine overall landscape forms and local relief details. Using vertical illumination, I see almost only the LRM colours, and by manipulating shaded relief illumination direction I can easily further enhance visibility of certain features. The colour-coding helps a lot to avoid errors due to optical illusions (i.e., apparent relief inversion).

Colour-coded LRM overlaid over a shaded relief image.

Colour-coded LRM overlaid over a shaded reliuef image. LIDAR data (c) LGL/LAD.

DEM visualisation techniques: Trend removal

In particular when using high-resolution DEM for archaeological prospection and mapping, we are interested in very subtle relief features: Many traces of former human activities (e.g. field boundaries) differ from their surroundings only by centimetres. Over periods of centuries and millenia, formerly distinct anthropogenic features such as burial mounds or fortifications become smoothed by erosion, infilling, ploughing or intentional levelling. As a result, archaeological relief features are often subtle in comparison with the natural topography on which they are superimposed.

Shaded relief image (illumination azimuth 315°, elevation 45°, no vertical exaggeration) showing a group of burial mounds, a rectangular Iron Age enclosure and modern field boundaries.

Shaded relief image (illumination azimuth 315°, elevation 45°, no vertical exaggeration) showing a group of burial mounds, a rectangular Iron Age enclosure and modern field boundaries. LIDAR data (c) LGL/LAD.

A simple method to extract small-scale relief features is trend removal. Just as in the two-dimensional case (e.g. removing the long-term trend from a time series to find anomalies), the three-dimensional case of removing the trend from a DEM is based on smoothing the data and then subtracting the smoothed data set from the original data set.

For smoothing the DEM, a (simple or weighted) average is computed for a moving window surrounding each pixel. For weighted averaging, a Gaussian function is commonly used. The intensity of DEM smoothing when applying such a low-pass filter depends on the size of the filter window and the standard deviation of the Gaussian filter. Subtracting the low-pass filtered DEM from the original DEM effectively means applying a high-pass filter.

The result is a raster map which contains positive and negative values for local bumps and depressions, respectively. To display this raster map as an image, greyscale or colour coding can be applied. Care has to be taken when interpreting such images. All convex and concave terrain edges are highlighted as relative positive and negative relief features, but this does not mean that they are positive and negative relief features in an absolute sense: the dark fringes surrounding the burial mounds in the example do not show actual ring ditches.

Greyscale image of a trend removal computed from the same DEM using a circular low pass with a radius of 12 pixels, values in the range from -1 to 1 are stretched to the black-white range). Note apparent ring ditches, in particular around the largest burial mound.

Greyscale image of a trend removal compouted from the same DEM using a circular low pass with a radius of 12 pixels, values in the range from -1 to 1 are stretched to the black-white range). Note apparent ring ditches, in particular around the largest burial mound. LIDAR data (c) LGL/LAD.

DEM visualisation techniques: Exaggerated relief

The directional illumination used in shaded relief visualisation can never show all relief details in an optimal way. Combining several shaded relief images computed for different illumination directions (e.g. by creating RGB colour composites or by applying Principal Component Analysis) can help, but the resulting images can be difficult to read.

Shaded relief image (illumination azimuth 315°, elevation 60°, no vertical exaggeration) of part of the Iron Age oppidum Heidengraben, Baden-Württemberg.

Shaded relief image (illumination azimuth 315°, elevation 60°, no vertical exaggeration) of part of the Iron Age oppidum Heidengraben, Baden-Württemberg. LIDAR data (c) LGL/LAD.

Because the optimal illumination for the recongnition of topographic detail depends to a large part on topographic position (mainly on slope and aspect), one option to improve relief feature visibility is to locally adapt illumination directions. Exaggerated Relief as proposed by Rusinkiewicz et al. (2006) is an algorithm that implements this idea. It’s not a simple process; it took me quite a while to fully understand how it works. Exaggerated Relief as described here is not to be confused with vertical exaggeration in creating shaded relief images.

It uses the same principle that underlies shaded relief illumination: pixel brightness is proportional to the cosine of the angle between the incoming light beam and the surface normal. But instead of using one global illumination direction for the whole DEM, illumination elevation is locally adjusted in a way that the light always strikes the surface at a shallow angle. Of course, doing this requires knowlwedge of surface slope and aspect; these parameters are computed for a smoothed (low-pass filtered) version of the same DEM. (To be precise, instead of smoothing the DEM, a surface normal map computed from the DEM is smoothed. This does not make a difference for the resulting images, but it makes computations easier.)

The resulting image for such a locally adapted illumination strongly emphasises local relief details, but it does not convey a visual impression of the overall topography. Therefore, the same computations are repeated at several scales: Another image for locally adapted illumination is calculated from the smoothed DEM, using a further, more strongly smoothed version of the same DEM to compute surface slope and aspect for locally adapting illumination elevation.

This process is repeated several times, smoothing the DEM surface normal map step by step and locally adjusting illumination elevation. Thereby, a series of images is created which each emphasise relief detail at a certain scale: the entire process is a multi-scale approach.

Pixel brightness for each of the resulting images is then multiplied by a constant exaggeration factor and limited (clamped) to a fixed range of values (-1…1). Finally, all images and one shaded relief image without local illumination adjustment (“base coat”) are combined to a single image by computing a weigthed average. To control how strongly the different images (i.e. the level of relief detail) influence the combined image, the weighting factor for each image is computed as a power of the standard deviation of the size of the Gaussian low pass filter which was used to smooth the DEM surface normal map and normalised so as to make the sum of weighting factoprs for all scales 1.

Several parameters can be adjusted to control the resulting image: global illumination direction (azimuth and elevation), minimum smoothing radius, number of scales, scale-to-scale incremental factor, exaggeration factor and weight exponent. To facilitate computation, a surface normal map (SNM) is computed from the DEM in a first processing step. Because all further computations are based on this SNM, other sources for SNMs (e.g. Polynomial Texture Mapping, PTM) can be used to compute Exaggerated Relief images.

Because Exaggerated Relief is based on the same fundamental principle of directional illumination, the resulting images are visually similar to shaded relief images. With appropriate settings, however, Exaggerated Relief conveys much more relief details largely irrespective of topographic position: subtle relief features are visible both in flat areas and on steep slopes, and there are almost no homogeneous white or black areas.

Exaggerated Relief visualisation of the same area (number of scales: 8, scale -to-scale factor: 1.414, exaggeration factor: 8, weight exponent: -0.5, global illumination azimuth: 315°, elevation: 45°). LIDAR data (c) LGL/LAD.

Exaggerated Relief visualisation of the same area (number of scales: 8, scale -to-scale factor: 1.414, exaggeration factor: 8, weight exponent: -0.5, global illumination azimuth: 315°, elevation: 45°). LIDAR data (c) LGL/LAD.

Exaggerated Relief visualisation of the same area (number of scales: 8, scale -to-scale factor: 1.414, exaggeration factor: 8, weight exponent: -0.5, global illumination azimuth: 315°, elevation: 45°). LIDAR data (c) LGL/LAD.

Exaggerated Relief visualisation of the same area (number of scales: 8, scale -to-scale factor: 1.414, exaggeration factor: 8, weight exponent: 0.0, global illumination azimuth: 315°, elevation: 45°). LIDAR data (c) LGL/LAD.

Exaggerated Relief visualisation of the same area (number of scales: 8, scale -to-scale factor: 1.414, exaggeration factor: 8, weight exponent: -0.5, global illumination azimuth: 315°, elevation: 45°). LIDAR data (c) LGL/LAD.

Exaggerated Relief visualisation of the same area (number of scales: 8, scale -to-scale factor: 1.414, exaggeration factor: 8, weight exponent: 0.5, global illumination azimuth: 315°, elevation: 45°). LIDAR data (c) LGL/LAD.

Does this mean that Exaggerated Relief is a perfect visualisation technique and that we don’t need any other technique? Not quite. One important issue is computation time. Because Exaggerated Relief is a relatively complicated algorithm working at multiple sucessive scales, computation times are high (increasing in particular with the number of scales).

Another problematic issue are visual artefacts resulting from the algorithm: depending on the settings (increasingly with a stronger influence of high-frequency relief details, i.e. a lower weight exponent), apparent banks are present along terrain edges. Interpretation of Exaggerated Relief images requires careful consideration of such artefacts and comparison with other visualisation techniques to avoid misinterpretations.

References

Rusinkiewicz, S., Burns, M., DeCarlo, D., 2006. Exaggerated Shading for depicting shape and detail. ACM Transactions on Graphics (Proceedings SIGGRAPH) 25(3), 1199–1205. [article link]