RoamesGeometry

Primitive geometric objects for Roames modelling.

Overview

This package provides a set of geometric primitives and containers for modelling objects in the physical world.

Points

The basic geometric point is 2D and 3D real, static vectors, like SVector{2, Float64} or SVector{3, Float32}. The convert2d and convert3d functions can be used to "upgrade" or "downgrade" dimensionality of points and other geometries.

There is functionality for loading and saving of pointclouds in either LAS or "Roames" HDF5 formats via load_pointcloud(filename) and save_pointcloud(filename, pointcloud).

Generally, in Julia v1.0 onwards a point cloud is represented as a Table from TypedTables.jl where the position column contains 3D points. For an introduction to using TypedTables see the user guide. Common manipulations on point clouds are shown below.

Geometries

RoamesGeometry supports a collection of 2D and 3D geometry types, common to many geospatial formats including ShapeFiles and Well-Known Text. Currently, we support:

• BoundingBox - a 3D region of space; an axis-aligned bounding box
• Line - a segment connecting two points
• LineString - a contiguous set of touching lines
• Quadratic - a quadratic function that supports "wind blown" wires
• Catenary - a vertically hanging catenary
• Polygon - a closed polygon (in the 2D sense), possibly with interior holes
• Sphere - a uniform sphere centered at a point with given radius
• Circle - a 2D circle
• Cylinder - A circle with min/max height
• Triangle - a three-sided polygon
• TriangularPrism - A triangle with min/max height

Generally, they can be transformed by a Transformation and can be plotted in Displaz via the plot3d function.

AbstractRegion and finding points

Some of these encompass 2D or 3D regions, and are subtypes of AbstractRegion. Regions are defined mathematically as open, connected and non-empty sets (or the closure thereof) and can support the in function (alternatitvely written in(point, region), point in region or point ∈ region).

For example, to ask the question if a point is inside a Sphere, one could write in(point, Sphere(centre, radius)). The in function supports "currying", so that f = in(Sphere(centre, radius)) is a function where f(point) returns true or false depending on whether point is in the Sphere or not. The predicate function f can then be used in higher order functions like map, filter and findall for filtering to just the points inside the Sphere.

To see this in action, given a pointcloud with a position column, we can construct a new pointcloud containing just the points inside a Sphere.

indices = findall(in(Sphere(centre, radius)), pointcloud.position)
pointcloud2 = pc[indices]

The findall function returns the indices of the elements matching a given predicate.

Lines and one-dimensional containers

Generally, points along one-dimensional primitives can be extracted via getindex, such as:

l = Line(SVector(0.0, 0.0, 0.0), SVector(2.0, 0.0, 0.0))
l[0.5] == SVector(0.5, 0.0, 0.0)

Note: geometries like Line can naturally be thought of as a collection of points, but they are not a region in the mathematic sense and calculation of point in line is not stable to e.g. floating-point rounding errors.

BoundingBox

Bounding box represents an axis-aligned rectangular prism and is an AbstractRegion. Its primary use is as a spatial acceleration structure, to check whether two objects lie within the same bounding box as an efficient pre-filtering step.

It supports the following interface:

• The boundingbox function is the primary constructor. You can call boundingbox with an arbitrary collection of geometries and a bounding box guaranteed to hold all the geometry is returned. Typically, it will be the smallest such bounding box. E.g. boundingbox(catenary, quadratic, line, point) or boundingbox(points::Vector{SVector{3,Float64}}).
• The in (or ∈) function/operator can indicate whether a point is within the bounding box, e.g. point ∈ box.
• The pad function extends the bounding box by a given amount, e.g. pad(bb, 1.0).
• The intersects function returns true if two bounding boxes interect, and false otherwise.
• The wireframe function returns the 12 lines outlining the box, as a Vector{Line{T}}.
• Displaz.plot3d can plot the wireframe of a bounding box directly.

Spatial acceleration using GridIndex

Spatial acceleration is used to speed up queries, like finding all the points within an AbstractRegion. A GridIndex is used to perform spatial acceleration on point cloud data of a "2.5 dimensional" nature - meaning points widely distributed in x and y, with a relatively few points in a given vertical column, such as those typical to aerial LiDAR.

A GridIndex tracks the points within each cell of an x-y grid of a given spacing, and furthermore orders them by height within the grid cell. This "index" can be used to make certain operations faster, for example to find all the points in a given region. Roughly speaking, given an AbstractRegion called region, spatial acceleration will

• Find bb = boundingbox(region).
• Use the grid to only search grid cells that intersect with bb.
• For each point in these cells, check precisely whether they are inside region.

That way, queries will skip the vast majority of points and the spatial index will make a tremendous performance improvement, often changing algorithms like PCA from O(n²) to O(n log n) or similar.

Acceleration indices are managed through the AcceleratedArrays.jl package. This package provides basic acceleration indices (like HashIndex and SortIndex) and is extended by RoamesGeometry to include GridIndex. An acceleration index is added to an array like so:

position = accelerate(position, GridIndex; spacing = 1.0)

Note that this has not mutated the original position array - rather it has created a new AcceleratedArray which wraps the old one. (Warning: mutating the positions will corrupt the index, meaning the results from findall and so-on will be incorrect).

For performance critical applications, one can re-order the array to be more cache-friendly and reduce lookups using the accelerate! function. If you do this to a point cloud, note the order of the other columns will not be modified, so the indices will get out of sync.

Distances

The distance between various geometries can be found. The distance function returns the (smallest) Euclidean distance between two geometric objects, and is currently defined between points and catenaries only. The closest_point and closest_points functions return the closest point(s) between geometries.

Distance to catenaries

The powerline_distances(catenary, point) function returns (dist, r, z, d), where:

• dist is the Euclidean distance to the closest point on the catenary.
• r is the horizontal distance perpendicular to the catenary.
• z is the height difference between the point and the catenary.
• d is the distance along the catenary, the end points being at catenary.lmin and catenary.lmax.

Input and Output

Well-known text

Input and output operations for well-known text is provided via the following functions:

• wkt(geometry) returns a String containing a well-known text representation of geometry.
• read_wkt(string) parses a WKT string and returns a geometry.
• load_wkt(filename) opens a well-known text file and reads a geometry.
• save_wkt(filename, geometry) saves geometry into a well-known text file.
• The lower-level operations read_wkt(io) and write_wkt(io, geometry) act on IO streams.

Point clouds

The load_pointcloud(filename) function can open .las and .h5 as a Table. To add a spatial acceleration GridIndex to the pointcloud, you must specify a grid spacing via load_pointcloud(filename, spacing = 1.0) (for a 1 metre grid). It is important to add the spatial index whenever you plan to make spatial queries, for example neighborhood search for noise filtering or PCA-based tasks.

The save_pointcloud(filename, pointcloud) function is able to save .las and .h5 files.

Both support GeoRepo2-style HDF5 files generated by ExtractPoints by specifying the relevant format string, e.g. load_pointcloud(filename, format = "XYZIrgb").

Example

To get started, consider this example where we combine the features of this library to classify points depending on whether they are inside a geometry or not, from files on disk.

using RoamesGeometry

pc = load_pointcloud("pointcloud.h5"; spacing = 1.0)
geom = load_wkt("geometry.wkt")

pc.classification .= 0
pc.classification[findall(in(geom)), pc.position)] .= 1

save_pointcloud("pointcloud2.h5", pc)