GitHub

Point Cloud Processing

Point cloud processing usually involves a number of steps which use the available point attributes to compute new attributes and eventually produce some derived data about the physical environment. PointClouds.jl includes the PointCloud struct for this purpose. It stores (a subset of) the original point attributes in a tabular format for in-memory processing and provides functionality to add/update/remove columns (i.e. point attributes) and to filter rows (i.e. points).

Loading point data

To load point data into memory, use the PointCloud constructor:

PointClouds.PointCloud — Type
PointCloud(input; kws...)

Load (a subset of) point data into memory for processing. The data is stored in “columnar format”, where each point attribute (coordinates, intensity, etc.) forms a column with one entry (row) per point. Attributes can be added/updated/removed in further processing steps.

The input can be one or multiple (passed as tuple/array) LAS values, or values that can be passed to read(input, LAS) such as the path to a LAS/LAZ file. Alternatively, a NamedTuple of Vectors can be passed, where the x, y, and z-names are interpreted as coordinates and all the other names as additional attributes.

Keywords

  • attributes: Additional attributes that are included for each point. Attributes are specified as a tuple/array of functions that are applied to each point to compute the attribute. By default, the name of the function is used as attribute name, but a manual name can be specified by passing a Pair{Symbol,Function} instead. The attributes can also be defined with a NamedTuple, where the keys are the attribute names and the values are the functions. Default: ().
  • coordinates: Subset of coordinates to load, e.g. (:x, :y) if the vertical coordinates are not required. Default: (:x, :y, :z).
  • crs: Coordinate reference system (CRS) that the x/y/z-coordinates should be transformed to, specified as any string understood by the PROJ library. Default: CRS of the first input.
  • x or lon: A tuple (xmin, xmax) with the minimum and maximum value of the x-coordinate range that should be retained, in the CRS of the output.
  • y or lat: A tuple (ymin, ymax) with the minimum and maximum value of the y-coordinate range that should be retained, in the CRS of the output.
  • z: A tuple (zmin, zmax) with the minimum and maximum value of the z-coordinate range that should be retained, in the CRS of the output.
  • filter: Predicate function that is called on each point to exclude points for which the filter function returns false.

Examples

PointCloud(las; attributes = (intensity, :t => gps_time))
PointCloud(las; attributes = (i = intensity, t = gps_time))
PointCloud(las; crs = "EPSG:4326", extent = ((-73.97, 40.80), (-73.95, 40.82)))
source

Coordinate reference system (CRS)

The getcrs function reads the CRS data from a PointCloud. Note that functions such as coordinates, filter, and the PointCloud constructor can also make use of this data without loading it explicitly.

Coordinates can be transformed when loading point data or with the transform function.

pts = transform(pts, crs = "EPSG:32618") # UTM 18N

This always produces a new variable, since the CRS is an immutable property of the PointCloud. Note that if the data previously had no CRS, the CRS is set to the target one but the coordinates are left untouched.

PointClouds.transform — Method
transform(p::PointCloud; crs)

Transform the coordinates of a PointCloud to a new coordinate reference system (CRS). If there is no current CRS, the coordinates remain unchanged but the new CRS is added to the PointCloud.

The required keyword argument crs can be any string understood by the PROJ library. If set to nothing, the CRS is removed.

source

Updating point attributes

Point attributes can be created, accessed, and updated using “property” syntax. The names of point attributes have no special meaning, except for x, y, and z, which are expected to represent the point coordinates. Some point processing functions may also attribute meaning to other names (e.g. neighbors!).

julia> pts = PointCloud((x = 0:2, y = 1:3));

julia> pts.x .+= 1
3-element Vector{Float64}:
 1.0
 2.0
 3.0

julia> pts.distance_from_origin = sqrt.(pts.x.^2 .+ pts.y.^2)
3-element Vector{Float64}:
 1.4142135623730951
 2.8284271247461903
 4.242640687119285

julia> keys(pts)
(:x, :y, :distance_from_origin)

julia> delete!(pts, :distance_from_origin);

Computing new attributes generally involves looping over the existing attributes. This can be done with the built-in map or broadcast functions, e.g. map((x, y) -> x^2 + y^2, pts.x, pts.y) or using the dot syntax as in the example above.

Parallel point processing

PointClouds.jl also includes the apply function, which is similar to map but runs in parallel using multi-threading and has some additional functionality specific to point-cloud processing.

PointClouds.Attributes.apply — Function
apply(f::Function, [T], p::PointCloud, attrs...; kws...)

Apply a function f to the attributes attrs of each point in p. The attribute names are passed as Symbols and the function should take one argument for each attribute name.

The element type of the output is determined automatically from the function and the argument types. It can also be set with the T argument in case the automatically determined type is too narrow or too wide.

Keywords

  • neighbors: Apply function to the attribute of a point’s neighbors rather than its own, if set to true or to an integer k. The function f receives AbstractVectors with the attribute values of the neighbors in that case. The indices of the k nearest neighbors are read from p.neighbors if available (see neighbors/neighbors!). If p.neighbors is unavailable or contains fewer than k indices, the neighbor search is (re)run (without updating p.neighbors). Default: false
source

Type stability

Warning

The PointCloud type contains no static information about the element type of attributes. Performance-sensitive point processing should therefore be done by a “kernel function” that takes individual attribute vectors as its argument. This is the case when using map, broadcast, or apply.

Coordinates can also be accessed with the coordinates function, which is type-stable:

PointClouds.coordinates — Method
coordinates(p::PointCloud, index; crs)

Obtain the x-, y-, and z-coordinate of the index-th point as a tuple of Float64s. The coordinate reference system (CRS) of the PointCloud is used unless a different crs is specified. To obtain coordinates of multiple points, pass a range of indices or : as the index argument.

Provide the crs keyword argument to transform the coordinates to a new CRS instead of the current CRS of the PointCloud.

source

Filtering points

The points in a PointCloud can be reduced to a subset with the filter/filter! functions. Alternatively, a subset can be selected through indexing (keepat! for in-place modification).

Base.filter — Method
filter([f::Function], pts::PointCloud, attrs...; kws...)

Create a copy of a PointCloud that excludes point which do not meet the specified filter criteria. Filtering can be done with a function f that returns false for each point that should be discarded, taking the attributes specified with attrs as arguments. Alternatively, points can be filtered by extent using the keyword arguments.

Note that you can also indexing (and keepat! for in-place modification) to filter a point-cloud, especially in combination with logical indexing.

Keywords

  • x or lon: A tuple (xmin, xmax) with the minimum and maximum value of the x-coordinate range that should be retained.
  • y or lat: A tuple (ymin, ymax) with the minimum and maximum value of the y-coordinate range that should be retained.
  • z: A tuple (zmin, zmax) with the minimum and maximum value of the z-coordinate range that should be retained.
  • crs: The coordinate reference system in which the extent should be applied. Set to the CRS of the PointCloud by default.

See also: filter!

source
Base.filter! — Method
filter!([f::Function], pts::PointCloud, attrs...; kws...)

Remove points that do not meet the specified filter criteria from a PointCloud. See filter for details.

source

Neighborhood-based processing

Several point-cloud processing algorithms depend on the neighborhood of each point. The neighbors function computes the indices of neighboring points, which can then be used in further processing steps, e.g. when using apply:

julia> pts = PointCloud((x = rand(16), y = rand(16), z = rand(16)));

julia> neighbors!(pts, 4);

julia> pts.zmin = apply(minimum, pts, :z; neighbors = true);

julia> pts.zavg = apply(zs -> sum(zs) / length(zs), pts, :z; neighbors = true);

Rasterization

Use rasterize to compute the mapping of points to a raster/voxel grid.

PointClouds.rasterize — Function
rasterize(p::PointCloud, dims; kws...)

Rasterize points in a PointCloud by assigning them to gridded locations within a target area. The raster has dimensions dims (given as a tuple of integers), is aligned with the axes of the x-/y-coordinates, and spans the full bounding box of the points, unless the extent is narrowed with keyword arguments. By default, each raster location is mapped to the points that lie within its “pixel area”, unless the radius or neighbors keyword arguments are specified. Currently, only 2D rasterization is supported.

Note that with the default “pixel area” method, each point gets assigned to exactly one raster location if it lies within the bounds of the raster. With the radius and neighbors methods, points may get assigned to multiple or zero raster locations. Furthermore, only the neighbors method guarantees that each raster location has points assigned to it; the other methods may produce empty pixels.

Keywords

  • x or lon: A tuple (xmin, xmax) with the minimum and maximum value of the x-coordinate range that should be divided into dims[1] intervals.
  • y or lat: A tuple (ymin, ymax) with the minimum and maximum value of the y-coordinate range that should be divided into dims[2] intervals.
  • radius: Assign all points within a distance of radius from the pixel center to that raster location. The unit of radius is the same as the x-/y-coordinates. Note that points may get assigned to multiple or zero raster locations.
  • neighbors: Assign a specific number neighbors of points to each raster location based on which ones are closest to the pixel center.
source

Mapping a function over the resulting RasterizedPointCloud applies it to each pixel/voxel, passing the list of points to the function.

zavg = map(r.z) do zs
    isempty(zs) ? missing : sum(zs) / length(zs)
end

This can also be used to combine raster and point data.

zstd = map(r.z, zavg) do zs, zavg
    length(zs) < 2 ? missing : sum(z -> abs2(z - zavg), zs) / (length(zs) - 1)
end