---
title: "Interfacing with `lidR`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Interfacing with `lidR`}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = requireNamespace("lidR")
)
```

The [`lidr` package](https://github.com/r-lidar/lidR) offers a robust suite of tools for processing LiDAR data. While `lacunr` does not require `lidR` as a strict dependency, it is assumed that most users will be working with point cloud data imported using `lidR`, and the package is designed to mesh well with `lidR`'s data objects. The following tips will help make combining these packages as seamless as possible.

```{r setup}
library(lacunr)
```

## Working with `LAS` objects

The workhorse of `lidR` is the `LAS` data object, an S4 class that emulates the structure of a .las/.laz point cloud file. `lacunr` does not include any example data in .las format, but we can generate a quick and dirty approximation by converting the built-in `glassfire` dataset to a `LAS` object:

```{r}
# convert 'glassfire' to LAS format
las_glassfire <- lidR::LAS(data = glassfire)
```

This triggers a message because we have not supplied all of the required metadata found in a proper .las file, but the new `LAS` object will suffice for the purposes of the following examples.

Now that `glassfire` has been converted, the point cloud `data.table` is now stored in an S4 slot called `@data`. If we try to feed the whole `LAS` object to `lacunr`'s `voxelize()` function, it will throw an error:

```{r error=TRUE}
# the wrong way to call voxelize() for a 'LAS' object:
vox <- voxelize(las_glassfire, edge_length = c(0.5, 0.5, 0.5))
```

Instead, we need to extract the `@data` slot and supply that to `voxelize()`:

```{r}
# voxelize the LAS point cloud, taking care to input the correct S4 slot
vox <- voxelize(las_glassfire@data, edge_length = c(0.5, 0.5, 0.5))
```

## Voxelization using `lidR`

`lidR` offers its own extremely versatile voxelization function, [`voxel_metrics()`](https://r-lidar.github.io/lidRbook/vba.html). This provides a useful alternative to `voxelize()`, although it is important to note that both functions utilize different algorithms and will not produce identical results (see the following section for more details). 

`voxel_metrics()` returns a `lasmetrics3d` object. `lacunr`'s `bounding_box()` function can accept this as an input, but it also requires that it contain a column named `N`, recording the number of points in each voxel. This column can be generated by `voxel_metrics()` using the following:

```{r}
# voxelize at 1m resolution, creating a column N containing the number of points
vox <- lidR::voxel_metrics(las_glassfire, ~list(N = length(Z)), res = 1)
# convert to array
box <- bounding_box(vox)
```

Failing to include this column will cause `bounding_box()` to throw an error.

## `lacunr::voxelize()` vs `lidR::voxel_metrics()`

`voxelize()` is adapted from the function `voxels()`, originally written by J. Antonio Guzmán Q. for the package [`rTLS`](https://github.com/Antguz/rTLS). It is intended as a complement rather than a replacement for `lidR`'s more elaborate `voxel_metrics()`. Each function has a different underlying algorithm and will produce distinct results from the same input data. The chief advantages of `voxelize()` over `voxel_metrics()` are: 

1. It allows --- and in fact requires --- the user to specify all three dimensions of the desired voxel resolution independently. This makes it possible to completely customize the shape of the voxels, in the rare instance that one wishes to divide up a point cloud into a non-cubic voxel grid. `voxel_metrics()` permits at most two dimensions.
2. The point cloud can be divided into an even number of voxel bins. For example, if you have a point cloud that spans 12 meters in the X dimension, and voxelize it at a resolution of 1 meter, the resulting data will be binned into 12 1-meter voxels along the X axis. The same point cloud will be binned into 13 voxels by `voxel_metrics()`. This is due to differences in how each function aligns the point cloud data within the voxel grid.