# 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
- Coordinate reference system (CRS)
- Updating point attributes
- Filtering points
- Neighborhood-based processing
- Rasterization

## 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 `Vector`

s 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)))
```

## 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.

`PointClouds.getcrs`

— Method`getcrs(pts::PointCloud)`

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.

## 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 `Symbol`

s 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`AbstractVector`

s 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`

### Type stability

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 `Float64`

s. 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`

.

## 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!`

`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.

## 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);
```

`PointClouds.Attributes.neighbors`

— Function`neighbors(p::PointCloud, k::Integer)`

Obtain the indices of the `k`

nearest neighbors of each point as a vector of vectors, where the inner vectors have a fixed length `k`

.

See also: `neighbors!`

`PointClouds.Attributes.neighbors!`

— Function`neighbors!(p::PointCloud, k::Integer)`

Store the indices of the `k`

nearest neighbors of each point in the `neighbors`

attribute of the `PointCloud`

.

See also: `neighbors`

## 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.

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
```