| Title: | Tools to Process Point Clouds Derived from Terrestrial Laser Scanning |
|---|---|
| Description: | A set of tools to process and calculate metrics on point clouds derived from terrestrial LiDAR (Light Detection and Ranging; TLS). Its creation is based on key aspects of the TLS application in forestry and ecology. Currently, the main routines are based on filtering, neighboring features of points, voxelization, canopy structure, and the creation of artificial stands. It is written using data.table and C++ language and in most of the functions it is possible to use parallel processing to speed-up the routines. |
| Authors: | J. Antonio Guzmán Q. [cre, aut, ths, cph] (ORCID: <https://orcid.org/0000-0002-0721-148X>), Ronny Hernandez [aut] (ORCID: <https://orcid.org/0000-0001-6225-7096>), Arturo Sanchez-Azofeifa [dgs, sad] (ORCID: <https://orcid.org/0000-0001-7768-6600>) |
| Maintainer: | J. Antonio Guzmán Q. <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.2.6.1 |
| Built: | 2026-05-20 09:49:37 UTC |
| Source: | https://github.com/antguz/rtls |
rTLS is a package that compiles a set of tools to process and calculate metrics on point clouds derived from terrestrial LiDAR (Light Detection and Ranging). Its creation is based on key aspects of the TLS application in forestry. Currently, the main routines are based on filtering, neighboring features of points, voxelization, optimal sphere or voxel size, and the creation of artificial stands. rTLS is written using data.table and C++ language and in most of the functions it is posible to use parallel processing to speed-up the routines.
Maintainer: J. Antonio Guzmán Q. [email protected] (ORCID) [thesis advisor, copyright holder]
Authors:
Ronny Hernandez [email protected] (ORCID)
Other contributors:
Arturo Sanchez-Azofeifa [email protected] (ORCID) [degree supervisor, scientific advisor]
Useful links:
Create an artificial forest stand of a given area using tree point clouds.
artificial_stand( files, n.trees, dimension, coordinates = NULL, sample = TRUE, replace = TRUE, overlap = NULL, rotation = TRUE, degrees = NULL, n_attempts = 100, progress = TRUE, plot = TRUE, ... )artificial_stand( files, n.trees, dimension, coordinates = NULL, sample = TRUE, replace = TRUE, overlap = NULL, rotation = TRUE, degrees = NULL, n_attempts = 100, progress = TRUE, plot = TRUE, ... )
files |
A |
n.trees |
A positive |
dimension |
A positive |
coordinates |
A |
sample |
Logical. If |
replace |
Logical. If |
overlap |
A positive |
rotation |
Logical. If |
degrees |
A positive |
n_attempts |
A positive |
progress |
Logical, if |
plot |
Logical. If |
... |
Parameters passed to |
When coordinates = NULL, artifical_stand adds, in sequence,
random coordinates to each files in the future stand based on the
crown area overlap. That is, first a tree from files is
randomly located within the stand dimention, then a second tree from
files will be located in the future stand based on the crown area
overlap from the previous tree, and so on. If during the random
location a given tree does not meet the requirements of overlap, new
random coordinates will be provided until the requirements are met.
Since artificial_stand will try to add tree to the stand until the
requirements are met, this could lead to an infinite loop if the stand
dimention is small or if the trees on files are large or many
n.trees. Therefore, the use of n_attempts is recommended to avoid
this scenario.
A list which contain a data.table (Trees) with the information of the point clouds used and their current coordinates in the stand, and another data.table with that compile all the point clouds used.
J. Antonio Guzmán Q.
#' #Import an example point cloud path <- system.file("extdata", "pc_tree.txt", package = "rTLS") #Creates a stand of 4 trees with 10% of overlap files <- rep(path, 4) artificial_stand(files, n.trees = 4, dimension = c(15, 15), overlap = 10) #Creates a stand of 4 trees with their locations location <- data.table(X = c(5, 10, 10, 5), Y = c(5, 5, 10, 10)) artificial_stand(files, n.trees = 4, dimension = c(15, 15), coordinates = location)#' #Import an example point cloud path <- system.file("extdata", "pc_tree.txt", package = "rTLS") #Creates a stand of 4 trees with 10% of overlap files <- rep(path, 4) artificial_stand(files, n.trees = 4, dimension = c(15, 15), overlap = 10) #Creates a stand of 4 trees with their locations location <- data.table(X = c(5, 10, 10, 5), Y = c(5, 5, 10, 10)) artificial_stand(files, n.trees = 4, dimension = c(15, 15), coordinates = location)
Estimates the canopy structure from a discrete returns scan from different TLS.
canopy_structure( TLS.type, scan, zenith.range, zenith.rings, azimuth.range, vertical.resolution, TLS.pulse.counts, TLS.resolution = NULL, TLS.coordinates = c(0, 0, 0), TLS.frame = NULL, TLS.angles = NULL, threads = 1 )canopy_structure( TLS.type, scan, zenith.range, zenith.rings, azimuth.range, vertical.resolution, TLS.pulse.counts, TLS.resolution = NULL, TLS.coordinates = c(0, 0, 0), TLS.frame = NULL, TLS.angles = NULL, threads = 1 )
TLS.type |
A |
scan |
If |
zenith.range |
If |
zenith.rings |
If |
azimuth.range |
A |
vertical.resolution |
A |
TLS.pulse.counts |
If |
TLS.resolution |
If |
TLS.coordinates |
A |
TLS.frame |
If |
TLS.angles |
A |
threads |
An |
Since scan describes discrete returns measured by the TLS, canopy_structre first simulates the number of pulses emitted based on Danson et al. (2007). The simulated pulses are
created based on the TLS properties (TLS.pulse.counts, TLS.resolution, TLS.frame) assuming that the scanner is perfectly balance. Then these pulses are rotated (rotate3D) based on the TLS.angles
roll, pitch, and yaw, and move to TLS.coordintates to simulate the positioning of the scanner during the scan. Rotated simulated-pulses of interest and scan returns are then extracted based on the zenith.range and azimuth.range for a given number of zenith.rings, azimuth.rings and vertical profiles.
The probability of gap (Pgap) is then estimated using the frequency of pulses and returns. For TLS.type = "multiple", the frequency of returns is estimated using the sum of 1/target count following Lovell et al. (2011).
Using the Pgap estimated per each zenith ring and vertical profile, canopy_structure then estimates the accumulative L(z) profiles based on the closest
zenith ring to 57.5 (hinge region) and, if TLS.type = "fixed.angle", the f(z) or commonly named PAVD based on the ratio of the
derivative of L(z) and height (z) following Jupp et al. 2009 (Equation 18). If TLS.type is equal to "single" or "multiple", canopy_structure also
estimates the normalized average weighted L/LAI, and then PAVD based on the L (hinge angle) at the highest height (LAI) and the ratio between the derivative
of L/LAI (average weighted) and the derivative of z (Jupp et al. 2009; Equation 21).
Jupp et al. 2009 excludes the zero zenith or fist ring to conduct the average weighted L/LAI estimations, canopy_structre does not excludes this sections since it depends on the regions of interest of the user.
Therefore, user should consider this difference since it may introduce more variability to profile estimations.
For any TLS.type, it returns a data.table with the height profiles defined by vertical.resolution, the gap probability based on the zenith.range and zenith.rings, and
the accumulative L(z) profiles based on the closest zenith ring to 57.5 degrees (hinge angle). If TLS.type is equal to "fixed.angle", it returns f(z) or commonly named PAVD based on
on the ratio of the derivative of L(z) and the derivative of height (z). If TLS.type is equal to "single" or "multiple", it returns the normalized average weighting L/LAI, and PAVD: based
on the L (hinge angle) at the highest height and the ratio between the derivative of L/LAI average weighted and the derivative of z.
J. Antonio Guzmán Q.
Danson F.M., Hetherington D., Morsdorf F., Koetz B., Allgower B. 2007. Forest canopy gap fraction from terrestrial laser scanning. IEEE Geosci. Remote Sensing Letters 4:157-160. doi: 10.1109/LGRS.2006.887064
Lovell J.L., Jupp D.L.B., van Gorsel E., Jimenez-Berni J., Hopkinson C., Chasmer L. 2011. Foliage profiles from ground based waveform and discrete point LiDAR. In: SilviLaser 2011, Hobart, Australia, 16–20 October 2011.
Jupp D.L.B., Culvenor D.S., Lovell J.L., Newnham G.J., Strahler A.H., Woodcock C.E. 2009. Estimating forest LAI profiles and structural parameters using a ground-based laser called “Echidna®”. Tree Physiology 29(2): 171-181. doi: 10.1093/treephys/tpn022
data(TLS_scan) #Using a multiple return file #Select the four columns required TLS_scan <- TLS_scan[, 1:4] #This will take a while# canopy_structure(TLS.type = "multiple", scan = TLS_scan, zenith.range = c(50, 70), zenith.rings = 4, azimuth.range = c(0, 360), vertical.resolution = 0.25, TLS.pulse.counts = c(2082, 580), TLS.frame = c(30, 130.024, 0, 359.90), TLS.angles = c(1.026, 0.760, -110.019)) #Using a single return file data(TLS_scan) #Subset to first return observations TLS_scan <- TLS_scan[Target_index == 1, 1:3] #This will take a while# canopy_structure(TLS.type = "single", scan = TLS_scan, zenith.range = c(50, 70), zenith.rings = 4, azimuth.range = c(0, 360), vertical.resolution = 0.25, TLS.pulse.counts = c(2082, 580), TLS.frame = c(30, 130.024, 0, 359.90), TLS.angles = c(1.026, 0.760, -110.019))data(TLS_scan) #Using a multiple return file #Select the four columns required TLS_scan <- TLS_scan[, 1:4] #This will take a while# canopy_structure(TLS.type = "multiple", scan = TLS_scan, zenith.range = c(50, 70), zenith.rings = 4, azimuth.range = c(0, 360), vertical.resolution = 0.25, TLS.pulse.counts = c(2082, 580), TLS.frame = c(30, 130.024, 0, 359.90), TLS.angles = c(1.026, 0.760, -110.019)) #Using a single return file data(TLS_scan) #Subset to first return observations TLS_scan <- TLS_scan[Target_index == 1, 1:3] #This will take a while# canopy_structure(TLS.type = "single", scan = TLS_scan, zenith.range = c(50, 70), zenith.rings = 4, azimuth.range = c(0, 360), vertical.resolution = 0.25, TLS.pulse.counts = c(2082, 580), TLS.frame = c(30, 130.024, 0, 359.90), TLS.angles = c(1.026, 0.760, -110.019))
Convert from East-North-Up cartesian coordinates to polar coordinates.
cartesian_to_polar(cartesian, anchor = c(0, 0, 0), digits = NULL)cartesian_to_polar(cartesian, anchor = c(0, 0, 0), digits = NULL)
cartesian |
A |
anchor |
A |
digits |
A |
It assumes that the positive *Z* axis is the reference vector for the zenith angle. Likewise, it assumes that the *Y* axis is the north-south direction (positive to negative) for the azimuth angle.
If a point from cartesian presents the same *XY* coordinates than anchor, angles returns NA.
A data.table with the zenith and azimuth angles (degrees), and the distance to the anchor coordinate.
J. Antonio Guzmán Q.
data(pc_tree) cartesian_to_polar(pc_tree) anchor <- c(1, 1, 1) cartesian_to_polar(pc_tree, anchor)data(pc_tree) cartesian_to_polar(pc_tree) anchor <- c(1, 1, 1) cartesian_to_polar(pc_tree, anchor)
Adaptive random sample consensus for cicle fitting.
circleRANSAC( cloud, fpoints, pconf, poutlier, max_iterations, threads = 1L, plot = TRUE )circleRANSAC( cloud, fpoints, pconf, poutlier, max_iterations, threads = 1L, plot = TRUE )
cloud |
A |
fpoints |
A |
pconf |
A |
poutlier |
A |
max_iterations |
An |
threads |
An |
plot |
Logical. If |
A data.table with the *XY* coordinate information of the circle center, the radius, the error based on the least squares fit, and the proportion of inliers.
J. Antonio Guzmán Q.
#Point cloud data("pc_tree") #Subset region at at breast height sub <- pc_tree[between(Z, 1.25, 1.35),] #Fit circle circleRANSAC(sub, fpoints = 0.2, pconf = 0.95, poutlier = c(0.5, 0.5), max_iterations = 100)#Point cloud data("pc_tree") #Subset region at at breast height sub <- pc_tree[between(Z, 1.25, 1.35),] #Fit circle circleRANSAC(sub, fpoints = 0.2, pconf = 0.95, poutlier = c(0.5, 0.5), max_iterations = 100)
Estimate the distance between a point and a group of point.
euclidean_distance(point, cloud, threads = 1L)euclidean_distance(point, cloud, threads = 1L)
point |
A |
cloud |
A |
threads |
An |
A numeric vector describing of point to each row of cloud.
J. Antonio Guzmán Q.
data("pc_tree") euclidean_distance(point = c(0, 0, 0), pc_tree)data("pc_tree") euclidean_distance(point = c(0, 0, 0), pc_tree)
Filtering of point clouds using different methods
filter( cloud, method, radius, min_neighbours, k, nSigma, edge_length, distance = "euclidean", threads = 1L, verbose = FALSE, progress = FALSE, ... )filter( cloud, method, radius, min_neighbours, k, nSigma, edge_length, distance = "euclidean", threads = 1L, verbose = FALSE, progress = FALSE, ... )
cloud |
A |
method |
A filtering method to use. It most be |
radius |
A |
min_neighbours |
An |
k |
An |
nSigma |
A |
edge_length |
A positive |
distance |
Type of distance to calculate. |
threads |
An |
verbose |
If TRUE, log messages to the console. |
progress |
If TRUE, log a progress bar when |
... |
Arguments passed to |
A data.table with the filtered points
J. Antonio Guzmán Q.
#Load data data("pc_tree") #Move pc_tree for comparison pc_compare <- pc_tree pc_compare$X <- pc_compare$X - 7 #SOR filter r1 <- filter(pc_tree, method = "SOR", k = 30, nSigma = 1) rgl::plot3d(r1, col = "red") #Filter rgl::points3d(pc_compare, col = "black") #Original #min_neighbours filter r2 <- filter(pc_tree, "min_neighbors", radius = 0.02, min_neighbours = 20) rgl::plot3d(r2, col = "red") #Filter rgl::points3d(pc_compare, col = "black") #Original #voxel_center filter r3 <- filter(pc_tree, method = "voxel_center", edge_length = 0.1) rgl::plot3d(r3, col = "red") #Filter rgl::points3d(pc_compare, col = "black") #Original#Load data data("pc_tree") #Move pc_tree for comparison pc_compare <- pc_tree pc_compare$X <- pc_compare$X - 7 #SOR filter r1 <- filter(pc_tree, method = "SOR", k = 30, nSigma = 1) rgl::plot3d(r1, col = "red") #Filter rgl::points3d(pc_compare, col = "black") #Original #min_neighbours filter r2 <- filter(pc_tree, "min_neighbors", radius = 0.02, min_neighbours = 20) rgl::plot3d(r2, col = "red") #Filter rgl::points3d(pc_compare, col = "black") #Original #voxel_center filter r3 <- filter(pc_tree, method = "voxel_center", edge_length = 0.1) rgl::plot3d(r3, col = "red") #Filter rgl::points3d(pc_compare, col = "black") #Original
Estimate geometry features of neighboring points in a cloud.
geometry_features( cloud, method, radius, k, max_neighbour, distance = "euclidean", target = FALSE, threads = 1L, verbose = FALSE, progress = TRUE, ... )geometry_features( cloud, method, radius, k, max_neighbour, distance = "euclidean", target = FALSE, threads = 1L, verbose = FALSE, progress = TRUE, ... )
cloud |
A |
method |
A character string specifying the method to estimate the neighbors. It most be one of |
radius |
A |
k |
An |
max_neighbour |
An |
distance |
Type of distance to calculate. |
target |
Logic. If |
threads |
An |
verbose |
If |
progress |
If |
... |
Arguments passed to |
The function returns the geometry features of the neighboring points
of a given point in cloud. Geometry features are represented by the
relative values of the eigenvalues derived from a covariance matrix of the
neighboring points. Geometry features are not estimated on target points
with less than 3 neighboring points.
A array describing the point of the cloud in rows,
the relative eigenvalues in columns, and the radius or k per slide.
If method = "radius_search", it add in the first column the number of
neighboring points.
J. Antonio Guzmán Q.
#Create cloud example <- data.table(X = runif(200, min=0, max=10), Y = runif(200, min=0, max=10), Z = runif(200, min=0, max=10)) #Using knn method with two different k k_test <- c(5, 10) geometry_features(example, method = "knn", k = k_test) #Using radius search method with two different radius radius_test <- c(3, 4) geometry_features(example, method = "radius_search", radius = radius_test, max_neighbour = 200)#Create cloud example <- data.table(X = runif(200, min=0, max=10), Y = runif(200, min=0, max=10), Z = runif(200, min=0, max=10)) #Using knn method with two different k k_test <- c(5, 10) geometry_features(example, method = "knn", k = k_test) #Using radius search method with two different radius radius_test <- c(3, 4) geometry_features(example, method = "radius_search", radius = radius_test, max_neighbour = 200)
Adapted K nearest neighbors based on RcppHNSW
knn( query, ref, k, distance = "euclidean", same = FALSE, threads = 1L, verbose = FALSE, progress = FALSE, ... )knn( query, ref, k, distance = "euclidean", same = FALSE, threads = 1L, verbose = FALSE, progress = FALSE, ... )
query |
A |
ref |
A |
k |
An |
distance |
Type of distance to calculate. |
same |
Logic. If |
threads |
An |
verbose |
If TRUE, log messages to the console. |
progress |
If TRUE, log a progress bar when |
... |
Arguments passed to |
This function is based on hnswlib C++ library (Malkov & Yashunin 2016) and
its bindings for R (RcppHNSW; Melville 2020) for a fast estimation of neighbors
points. It is adapted to simplify the workflow within rTLS.
If you use this function, please consider cite the C++ library and
RcppHNSW package.
A data.table with three columns describing the indices of the query, ref, and k neighbors and the distances.
J. Antonio Guzmán Q.
Malkov, Y. A., & Yashunin, D. A. (2016). Efficient and robust approximate nearest neighbor search using Hierarchical Navigable Small World graphs. arXiv preprint arXiv:1603.09320.
#Point cloud data("pc_tree") #knn search using k = 3 knn(pc_tree, pc_tree, k = 3, same = TRUE)#Point cloud data("pc_tree") #knn search using k = 3 knn(pc_tree, pc_tree, k = 3, same = TRUE)
Intersection of a line by an Axis-Aligned Bounding Box.
line_AABB(orig, end, AABB_min, AABB_max)line_AABB(orig, end, AABB_min, AABB_max)
orig |
A |
end |
A |
AABB_min |
A |
AABB_max |
A |
The interaction of a line with a AABB may result in five scenarios:
i) the line is not intercepted by a AAABB (0), ii) the origin and end
of the line falls within the AABB (1), iii) the origin point of the
line falls within the AABB both not the end point (2), iv) the end point
of the line falls within the AABB both not the origin point (3), and
v) the line is intercepted by the AABB (4).
An numeric vector of length two, describing if the line was
intercepted or not, and the length of the intercepted line within in the AABB.
See details.
J. Antonio Guzmán Q.
#Create origins and end paths orig <- data.table(X = c(0, 0, 0, 0, 0), Y = c(-0.45, -0.25, 0, 0.25, 0.45), Z = c(-1, -0.25, 0, -1, -1)) end <- data.table(X = c(0, 0, 0, 0, 0), Y = c(-0.45, -0.25, 0, 0.25, 0.45), Z = c(-0.75, 0.25, 1, 0, 1)) #Create the AABB AABB <- matrix(c(0, 0, 0), ncol = 3) edge_length <- c(1, 1, 1) AABB_min <- c(AABB[1, 1] - edge_length[1]/2, AABB[1, 2] - edge_length[2]/2, AABB[1, 3] - edge_length[3]/2) AABB_max <- c(AABB[1, 1] + edge_length[1]/2, AABB[1, 2] + edge_length[2]/2, AABB[1, 3] + edge_length[3]/2) #Plot cube <- rgl::cube3d() cube <- rgl::scale3d(cube, edge_length[1]/2, edge_length[2]/2, edge_length[3]/2) box <- rgl::translate3d(cube, AABB[1, 1], AABB[1, 2], AABB[1, 3]) rgl::shade3d(box, col= "green", alpha = 0.6) rgl::points3d(orig, size = 4, col = "black") rgl::points3d(end, size = 4, col = "red") #Line no intercepted rgl::lines3d(c(orig[1, 1], end[1, 1]), c(orig[1, 2], end[1, 2]), c(orig[1, 3], end[1, 3]), col = "grey") line_AABB(orig[1,], end[1,], AABB_min, AABB_max) #Both ends falls inside rgl::lines3d(c(orig[2, 1], end[2, 1]), c(orig[2, 2], end[2, 2]), c(orig[2, 3], end[2, 3]), col = "red") line_AABB(orig[2,], end[2,], AABB_min, AABB_max) #Oring falls inside, but not the end. rgl::lines3d(c(orig[3, 1], end[3, 1]), c(orig[3, 2], end[3, 2]), c(orig[3, 3], end[3, 3]), col = "blue") line_AABB(orig[3,], end[3,], AABB_min, AABB_max) #End falls inside, but not the orig rgl::lines3d(c(orig[4, 1], end[4, 1]), c(orig[4, 2], end[4, 2]), c(orig[4, 3], end[4, 3]), col = "green") line_AABB(orig[4,], end[4,], AABB_min, AABB_max) #Some segments of the line are intercepted rgl::lines3d(c(orig[5, 1], end[5, 1]), c(orig[5, 2], end[5, 2]), c(orig[5, 3], end[5, 3]), col = "black") line_AABB(orig[5,], end[5,], AABB_min, AABB_max)#Create origins and end paths orig <- data.table(X = c(0, 0, 0, 0, 0), Y = c(-0.45, -0.25, 0, 0.25, 0.45), Z = c(-1, -0.25, 0, -1, -1)) end <- data.table(X = c(0, 0, 0, 0, 0), Y = c(-0.45, -0.25, 0, 0.25, 0.45), Z = c(-0.75, 0.25, 1, 0, 1)) #Create the AABB AABB <- matrix(c(0, 0, 0), ncol = 3) edge_length <- c(1, 1, 1) AABB_min <- c(AABB[1, 1] - edge_length[1]/2, AABB[1, 2] - edge_length[2]/2, AABB[1, 3] - edge_length[3]/2) AABB_max <- c(AABB[1, 1] + edge_length[1]/2, AABB[1, 2] + edge_length[2]/2, AABB[1, 3] + edge_length[3]/2) #Plot cube <- rgl::cube3d() cube <- rgl::scale3d(cube, edge_length[1]/2, edge_length[2]/2, edge_length[3]/2) box <- rgl::translate3d(cube, AABB[1, 1], AABB[1, 2], AABB[1, 3]) rgl::shade3d(box, col= "green", alpha = 0.6) rgl::points3d(orig, size = 4, col = "black") rgl::points3d(end, size = 4, col = "red") #Line no intercepted rgl::lines3d(c(orig[1, 1], end[1, 1]), c(orig[1, 2], end[1, 2]), c(orig[1, 3], end[1, 3]), col = "grey") line_AABB(orig[1,], end[1,], AABB_min, AABB_max) #Both ends falls inside rgl::lines3d(c(orig[2, 1], end[2, 1]), c(orig[2, 2], end[2, 2]), c(orig[2, 3], end[2, 3]), col = "red") line_AABB(orig[2,], end[2,], AABB_min, AABB_max) #Oring falls inside, but not the end. rgl::lines3d(c(orig[3, 1], end[3, 1]), c(orig[3, 2], end[3, 2]), c(orig[3, 3], end[3, 3]), col = "blue") line_AABB(orig[3,], end[3,], AABB_min, AABB_max) #End falls inside, but not the orig rgl::lines3d(c(orig[4, 1], end[4, 1]), c(orig[4, 2], end[4, 2]), c(orig[4, 3], end[4, 3]), col = "green") line_AABB(orig[4,], end[4,], AABB_min, AABB_max) #Some segments of the line are intercepted rgl::lines3d(c(orig[5, 1], end[5, 1]), c(orig[5, 2], end[5, 2]), c(orig[5, 3], end[5, 3]), col = "black") line_AABB(orig[5,], end[5,], AABB_min, AABB_max)
Intersection of lines by several Axis-Aligned Bounding Boxs.
lines_interception(orig, end, AABBs, edge_length, threads = 1, progress = TRUE)lines_interception(orig, end, AABBs, edge_length, threads = 1, progress = TRUE)
orig |
A |
end |
A |
AABBs |
A |
edge_length |
A positive |
threads |
An |
progress |
Logical, if |
It returns a data.table with nine columns: 1-5 columns with
the counts for the code of intersection (see line_AABB), and
6-9 columns with sum the path length of intersection. The number of rows match
with nrow(AABBs).
J. Antonio Guzmán Q.
#Create points with paths n <- 20 orig <- data.table(X = runif(n, min = -5, max = 5), Y = runif(n, min = -5, max = 5), Z = runif(n, min = -5, max = 5)) end <- data.table(X = runif(n, min = -5, max = 5), Y = runif(n, min = -5, max = 5), Z = runif(n, min = -5, max = 5)) #Create a potential AABB AABBs <- data.table(X = 0, Y = 0, Z = 0) edge_length <- c(2, 2, 2) #Plot cube <- rgl::cube3d() cube <- rgl::scale3d(cube, edge_length[1]/2, edge_length[2]/2, edge_length[3]/2) box <- rgl::translate3d(cube, AABBs[[1]], AABBs[[2]], AABBs[[3]]) rgl::shade3d(box, col= "green", alpha = 0.6) rgl::points3d(orig, size = 5, col = "black") rgl::points3d(end, size = 5, col = "red") for(i in 1:nrow(orig)) { rgl::lines3d(c(orig[[1]][i], end[[1]][i]), c(orig[[2]][i], end[[2]][i]), c(orig[[3]][i], end[[3]][i]), col = "grey") } #Estimation lines_interception(orig, end, AABBs, edge_length, progress = FALSE)#Create points with paths n <- 20 orig <- data.table(X = runif(n, min = -5, max = 5), Y = runif(n, min = -5, max = 5), Z = runif(n, min = -5, max = 5)) end <- data.table(X = runif(n, min = -5, max = 5), Y = runif(n, min = -5, max = 5), Z = runif(n, min = -5, max = 5)) #Create a potential AABB AABBs <- data.table(X = 0, Y = 0, Z = 0) edge_length <- c(2, 2, 2) #Plot cube <- rgl::cube3d() cube <- rgl::scale3d(cube, edge_length[1]/2, edge_length[2]/2, edge_length[3]/2) box <- rgl::translate3d(cube, AABBs[[1]], AABBs[[2]], AABBs[[3]]) rgl::shade3d(box, col= "green", alpha = 0.6) rgl::points3d(orig, size = 5, col = "black") rgl::points3d(end, size = 5, col = "red") for(i in 1:nrow(orig)) { rgl::lines3d(c(orig[[1]][i], end[[1]][i]), c(orig[[2]][i], end[[2]][i]), c(orig[[3]][i], end[[3]][i]), col = "grey") } #Estimation lines_interception(orig, end, AABBs, edge_length, progress = FALSE)
Estimate the minimum distance between points in a point cloud.
min_distance( cloud, distance = "euclidean", threads = 1L, verbose = FALSE, progress = FALSE, ... )min_distance( cloud, distance = "euclidean", threads = 1L, verbose = FALSE, progress = FALSE, ... )
cloud |
A |
distance |
Type of distance to calculate. |
threads |
An |
verbose |
If TRUE, log messages to the console. |
progress |
If TRUE, log a progress bar when |
... |
Arguments passed to |
A numeric vector describing the minimum distance between points.
J. Antonio Guzmán Q.
data("pc_tree") #Estimate the minimum distance of a sample o 100 points min_distance(pc_tree)data("pc_tree") #Estimate the minimum distance of a sample o 100 points min_distance(pc_tree)
A data.table from a point cloud of a tree with a spatial point resolution of 0.05 mm.
data(pc_tree)data(pc_tree)
A data.table with three columns, which are:
the "X" coordinate
the "Y coordinate
the "Z" coordinate
A data.table where the rows represent the points and the three columns represent the *XYZ* coordinates.
Guzman, Sharp, Alencastro, Sanchez-Azofeifa. 2018. To be published.
data(pc_tree) head(pc_tree)data(pc_tree) head(pc_tree)
The plot method for objects of class "voxels" created using the voxels function.
plot_voxels( voxels, add.points = TRUE, add.voxels = TRUE, border = TRUE, points.size = 1, points.col = "black", fill.col = "forestgreen", line.lwd = 0.5, line.col = "black", alpha = 0.1, ... )plot_voxels( voxels, add.points = TRUE, add.voxels = TRUE, border = TRUE, points.size = 1, points.col = "black", fill.col = "forestgreen", line.lwd = 0.5, line.col = "black", alpha = 0.1, ... )
voxels |
Object of class |
add.points |
Logical, if |
add.voxels |
Logical, if |
border |
Logical, if |
points.size |
The points size, a positive number to use if plot |
points.col |
A |
fill.col |
A |
line.lwd |
The line width, a positive number, defaulting to 0.5. |
line.col |
A |
alpha |
A positive numeric vector describing the transparency of the voxels to fill. This value most be between 0.0 (fully transparent) .. 1.0 (opaque). |
... |
General arguments passed to |
A 3D plot of a point cloud and voxels.
J. Antonio Guzmán Q.
voxels, voxels_counting, summary_voxels
data("pc_tree") ###Create cubes of a size of 7x7x3.5. vox <- voxels(pc_tree, edge_length = c(7, 7, 3.5)) plot_voxels(vox)data("pc_tree") ###Create cubes of a size of 7x7x3.5. vox <- voxels(pc_tree, edge_length = c(7, 7, 3.5)) plot_voxels(vox)
Convert from polar to cartesian coordinates.
polar_to_cartesian(polar, threads = 1, digits = NULL)polar_to_cartesian(polar, threads = 1, digits = NULL)
polar |
A |
threads |
An |
digits |
A |
A data.table with three columns describing the *XYZ* of the cartesian coordinates.
J. Antonio Guzmán Q.
#Creates a hemisphere of points each 2 degrees zenith <- seq(0, 90, 2) azimuth <- seq(0, 360, 2) hemi <- CJ(zenith, azimuth) hemi$distance <- 1 hemicloud <- polar_to_cartesian(hemi) rgl::plot3d(hemicloud)#Creates a hemisphere of points each 2 degrees zenith <- seq(0, 90, 2) azimuth <- seq(0, 360, 2) hemi <- CJ(zenith, azimuth) hemi$distance <- 1 hemicloud <- polar_to_cartesian(hemi) rgl::plot3d(hemicloud)
Adapted radius searching of points based on RcppHNSW
radius_search( query, ref, radius, max_neighbour, distance = "euclidean", same = FALSE, threads = 1L, verbose = FALSE, progress = FALSE, ... )radius_search( query, ref, radius, max_neighbour, distance = "euclidean", same = FALSE, threads = 1L, verbose = FALSE, progress = FALSE, ... )
query |
A |
ref |
A |
radius |
A |
max_neighbour |
An |
distance |
Type of distance to calculate. |
same |
Logic. If |
threads |
An |
verbose |
If TRUE, log messages to the console. |
progress |
If TRUE, log a progress bar when |
... |
Arguments passed to |
This function is based on hnswlib C++ library (Malkov & Yashunin 2016) and
its bindings for R (RcppHNSW; Melville 2020) for a fast estimation of neighbors
points. It is adapted to simplify the workflow within rTLS.
If you use this function, please consider cite the C++ library and
RcppHNSW package.
A data.table with three columns describing the indices of the query and ref points and the distances.
J. Antonio Guzmán Q.
Malkov, Y. A., & Yashunin, D. A. (2016). Efficient and robust approximate nearest neighbor search using Hierarchical Navigable Small World graphs. arXiv preprint arXiv:1603.09320.
#Point cloud data("pc_tree") #Radius search of 1 radius_search(pc_tree, pc_tree, radius = 1, max_neighbour = 100)#Point cloud data("pc_tree") #Radius search of 1 radius_search(pc_tree, pc_tree, radius = 1, max_neighbour = 100)
Rotate a plane of coordinates to a given angle.
rotate2D(plane, angle, threads = 1)rotate2D(plane, angle, threads = 1)
plane |
A |
angle |
A |
threads |
An |
A data.table with the rotation applied to plane.
J. Antonio Guzmán Q.
data(pc_tree) plot(pc_tree[,1:2]) #Rotate in 45 degrees using Z axis of the cloud plot(rotate2D(pc_tree[,1:2], angle = 45))data(pc_tree) plot(pc_tree[,1:2]) #Rotate in 45 degrees using Z axis of the cloud plot(rotate2D(pc_tree[,1:2], angle = 45))
Rotate point cloud based on the roll, pitch, and yaw angles.
rotate3D(cloud, roll = 0, pitch = 0, yaw = 0, threads = 1)rotate3D(cloud, roll = 0, pitch = 0, yaw = 0, threads = 1)
cloud |
A |
roll |
A |
pitch |
A |
yaw |
A |
threads |
An |
The *XYZ* coordinates are transformed to E-N-U coordinates (ENU system, East-North-Up).
A data.table with the rotation applied to cloud.
J. Antonio Guzmán Q.
data(pc_tree) rgl::plot3d(pc_tree) rgl::plot3d(rotate3D(pc_tree, roll = 45, pitch = 45, yaw = 0))data(pc_tree) rgl::plot3d(pc_tree) rgl::plot3d(rotate3D(pc_tree, roll = 45, pitch = 45, yaw = 0))
Applies the voxels_counting function on a grid base point cloud.
stand_counting( cloud, xy.res, z.res = NULL, points.min = NULL, min_size, edge_sizes = NULL, length_out = 10, bootstrap = FALSE, R = NULL, progress = TRUE, parallel = FALSE, threads = NULL )stand_counting( cloud, xy.res, z.res = NULL, points.min = NULL, min_size, edge_sizes = NULL, length_out = 10, bootstrap = FALSE, R = NULL, progress = TRUE, parallel = FALSE, threads = NULL )
cloud |
A |
xy.res |
A positive |
z.res |
A positive |
points.min |
A positive |
min_size |
A positive |
edge_sizes |
A positive |
length_out |
A positive |
bootstrap |
Logical. If |
R |
A positive |
progress |
Logical, if |
parallel |
Logical, if |
threads |
An |
A data.table with the summary of the voxels per grid created with their features.
J. Antonio Guzmán Q.
voxels_counting, voxels, summary_voxels
data(pc_tree) #Applying stand_counting. stand_counting(pc_tree, xy.res = c(4, 4), min_size = 3) #Applying stand_counting using bootstrap in the H index. stand_counting(pc_tree, xy.res = c(4, 4), min_size = 3, bootstrap = TRUE, R = 10)data(pc_tree) #Applying stand_counting. stand_counting(pc_tree, xy.res = c(4, 4), min_size = 3) #Applying stand_counting using bootstrap in the H index. stand_counting(pc_tree, xy.res = c(4, 4), min_size = 3, bootstrap = TRUE, R = 10)
Create a summary objects of class "voxels" created using the voxels.
summary_voxels(voxels, edge_length = NULL, bootstrap = FALSE, R = NULL)summary_voxels(voxels, edge_length = NULL, bootstrap = FALSE, R = NULL)
voxels |
An object of class |
edge_length |
A positive |
bootstrap |
Logical, if |
R |
A positive |
The function provides 12 main statistics of the voxels. Specifically, the first three columns represent the edge length of the voxels, the following three columns (ei. N_voxels, Volume, Surface) describe the number of voxels created, the total volume that they represent, and the surface area that they cover.
Following columns represent the mean (Density_mean) and sd (Density_sd) of the density of points per voxel (e.g. points/m2). Columns 9:12 provide metrics calculated using the Shannon Index. Specifically, H describe the entropy, H_max the maximum entropy, Equitavility the ratio between H and Hmax, and Negentropy describe the product of Hmax - H.
If bootstrap = TRUE four more columns are created (13:16). These represent the mean and sd of the H index estimated using bootstrap (H_boot_mean and H_boot_sd), the Equtavility_boot as the ratio of the ratio between H_boot_sd and Hmax, and Negentropy_boot as the product Hmax - H_boot_mean.
A data.table with with the summary of voxels.
J. Antonio Guzmán Q.
voxels, voxels_counting, plot_voxels
data("pc_tree") #Apply a summary on a object of class "voxels" using bootstrap with 1000 replicates. vox <- voxels(pc_tree, edge_length = c(0.5, 0.5, 0.5)) summary_voxels(vox, bootstrap = TRUE, R = 1000) #Apply a summary on a product from 'voxels' using bootstrap with 1000 replicates. vox <- voxels(pc_tree, edge_length = c(0.5, 0.5, 0.5), obj.voxels = FALSE) summary_voxels(vox, edge_length = c(0.5, 0.5, 0.5), bootstrap = TRUE, R = 1000)data("pc_tree") #Apply a summary on a object of class "voxels" using bootstrap with 1000 replicates. vox <- voxels(pc_tree, edge_length = c(0.5, 0.5, 0.5)) summary_voxels(vox, bootstrap = TRUE, R = 1000) #Apply a summary on a product from 'voxels' using bootstrap with 1000 replicates. vox <- voxels(pc_tree, edge_length = c(0.5, 0.5, 0.5), obj.voxels = FALSE) summary_voxels(vox, edge_length = c(0.5, 0.5, 0.5), bootstrap = TRUE, R = 1000)
A data.table from a TLS scan.
data(TLS_scan)data(TLS_scan)
A data.table with five columns, which are:
the "X" coordinate
the "Y coordinate
the "Z" coordinate
The number of received by the same laser shot
The rank of the returned pulse in the target count of received by the same laser shot
A data.table where the rows represent the pulse returns and the
three columns represent the *XYZ* coordinates, and the target count and index.
A TLS scan conducted using a Reigel VZ400i with a vertical and
horizontal resolution of 0.048 and 0.622 degrees (2082 and 580 lines, respectively).
The scanner has frame of scanning between 30 and 130.024 degrees zenith and 0 and
359.90 degrees azimuth. At the moment of the scan the roll, pitch, and yaw of
the scanner were 1.026, 0.746, -110.019, respectively. The scanner coordinates
in this scan are x = 0, y = 0, z = 0.
data(TLS_scan) head(TLS_scan)data(TLS_scan) head(TLS_scan)
Estimate the tree height, crown area, and the diameter at breast height of a tree point cloud
tree_metrics(cloud, region.diameter = NULL, relocateZ = TRUE)tree_metrics(cloud, region.diameter = NULL, relocateZ = TRUE)
cloud |
A |
region.diameter |
A |
relocateZ |
Logical, if |
The tree height is estimated based on the maximum value of *Z*, the
crown area is calculated applying a convex hull on the point cloud, while the
DBH is calculated extracting the area of the convex hull on the subset of points
between region.diameter, and then estimating the diameter of a circle.
For another estimation of DBH try circleRANSAC or for irregular
trucks try trunk_volume.
A data.table with the tree height, crown area, and diameter
J. Antonio Guzman Q. and Ronny Hernandez
data("pc_tree") tree_metrics(pc_tree)data("pc_tree") tree_metrics(pc_tree)
Estimates the tree trunk volume of a point cloud using the ashape3d package.
trunk_volume(cloud, max.height = NULL, alpha = 0.2, plot = TRUE, ...)trunk_volume(cloud, max.height = NULL, alpha = 0.2, plot = TRUE, ...)
cloud |
A |
max.height |
A |
alpha |
A |
plot |
Logical. If |
... |
General arguments passed to |
This is an adaptation of the code develop by Lafarge & Pateiro-Lopez (2017) based on Edelsbrunner & Mucke (1994) for the quick extraction of the tree trunk volume. Therefore, if you use this code we kindly suggest to cite these documents in your research.
A numeric vector with the estimated trunk volume.
J. Antonio Guzmán Q.
Lafarge, T., Pateiro-Lopez, B. (2017). Implementation of the 3D Alpha-Shape for the Reconstruction of 3D Sets from a Point Cloud. Available at https://CRAN.R-project.org/package=alphashape3d.
Edelsbrunner, H., Mucke, E. P. (1994). Three-Dimensional Alpha Shapes. ACM Transactions on Graphics, 13(1), pp.43-72.
data("pc_tree") #Estimates the trunk volume of a height lower than 1.75. trunk_volume(pc_tree, max.height = 1.75)data("pc_tree") #Estimates the trunk volume of a height lower than 1.75. trunk_volume(pc_tree, max.height = 1.75)
Create cubes of a given distance in a point cloud though their voxelization. It use a modify version of the code used in Greaves et al. 2015.
voxels(cloud, edge_length, threads = 1L, obj.voxels = TRUE)voxels(cloud, edge_length, threads = 1L, obj.voxels = TRUE)
cloud |
A |
edge_length |
A positive |
threads |
An |
obj.voxels |
Logical. If |
Voxels are created from the negative to the positive *XYZ* coordinates.
If obj.voxels == TRUE, it return an object of class "voxels" which contain a list with the points used to create the voxels, the parameter edge_length, and the voxels created. If FALSE, it returns a data.table with the coordinates of the voxels created and the number of points in each voxel.
J. Antonio Guzmán Q.
Greaves, H. E., Vierling, L. A., Eitel, J. U., Boelman, N. T., Magney, T. S., Prager, C. M., & Griffin, K. L. (2015). Estimating aboveground biomass and leaf area of low-stature Arctic shrubs with terrestrial LiDAR. Remote Sensing of Environment, 164, 26-35.
voxels_counting, plot_voxels, summary_voxels
data("pc_tree") ###Create cube of a size of 0.5. voxels(pc_tree, edge_length = c(0.5, 0.5, 0.5))data("pc_tree") ###Create cube of a size of 0.5. voxels(pc_tree, edge_length = c(0.5, 0.5, 0.5))
Creates cube like voxels of different size on a point cloud using the voxels function, and then return a summary_voxels of their features.
voxels_counting( cloud, edge_sizes = NULL, min_size, length_out = 10, bootstrap = FALSE, R = NULL, progress = TRUE, parallel = FALSE, threads = NULL )voxels_counting( cloud, edge_sizes = NULL, min_size, length_out = 10, bootstrap = FALSE, R = NULL, progress = TRUE, parallel = FALSE, threads = NULL )
cloud |
A |
edge_sizes |
A positive |
min_size |
A positive |
length_out |
A positive |
bootstrap |
Logical. If |
R |
A positive |
progress |
Logical, if |
parallel |
Logical, if |
threads |
An |
A data.table with the summary of the voxels created with their features.
J. Antonio Guzmán Q.
voxels, summary_voxels, plot_voxels
data(pc_tree) #Applying voxels counting. voxels_counting(pc_tree, min_size = 2) #Voxels counting using bootstrap on the H indexes with 1000 repetitions. voxels_counting(pc_tree, min_size = 2, bootstrap = TRUE, R = 1000)data(pc_tree) #Applying voxels counting. voxels_counting(pc_tree, min_size = 2) #Voxels counting using bootstrap on the H indexes with 1000 repetitions. voxels_counting(pc_tree, min_size = 2, bootstrap = TRUE, R = 1000)