Nothing
R/stem_points_methods.R
In TreeLS: Terrestrial Point Cloud Processing of Forest Data
Defines functions
stm.eigen.knn
stm.hough
Documented in
stm.eigen.knn
stm.hough
# ===============================================================================
#
# Developers:
#
# Tiago de Conto - tdc.florestal@gmail.com - https://github.com/tiagodc/
#
# COPYRIGHT: Tiago de Conto, 2020
#
# This piece of software is open and free to use, redistribution and modifications
# should be done in accordance to the GNU General Public License >= 3
# Use this software as you wish, but no warranty is provided whatsoever. For any
# comments or questions on TreeLS, please contact the developer (prefereably through my github account)
#
# If publishing any work/study/research that used the tools in TreeLS,
# please don't forget to cite the proper sources!
#
# Enjoy!
#
# ===============================================================================
#' Stem denoising algorithm: Hough Transform
#' @description This function is meant to be used inside \code{\link{stemPoints}}. It applies an adapted version of the Hough Transform for circle search. Mode details are given in the sections below.
#' @template param-h_step
#' @template param-max-d
#' @template param-hbase
#' @template param-pixel-size
#' @template param-min-density
#' @template param-min-votes
#' @section \code{LAS@data} Special Fields:
#'
#' Meaninful new fields in the output:
#'
#' \itemize{
#' \item \code{Stem}: \code{TRUE} for stem points
#' \item \code{Segment}: stem segment number (from bottom to top and nested with TreeID)
#' \item \code{Radius}: approximate radius of the point's stem segment estimated by the Hough Transform - always a multiple of the \code{pixel_size}
#' \item \code{Votes}: votes received by the stem segment's center through the Hough Transform
#' }
#'
#' @template section-hough-transform
#' @template reference-olofsson
#' @template reference-thesis
#' @export
stm.hough = function(h_step=0.5, max_d=0.5, h_base = c(1,2.5), pixel_size=0.025, min_density=0.1, min_votes=3){
if(length(h_base) != 2)
stop('h_base must be a numeric vector of length 2')
if(diff(h_base) <= 0)
stop('h_base[2] must be larger than h_base[1]')
params = list(
h_step = h_step,
max_d = max_d,
pixel_size = pixel_size,
min_density = min_density,
min_votes = min_votes
)
for(i in names(params)){
val = params[[i]]
if(length(val) != 1)
stop( i %>% paste('must be of length 1') )
if(!is.numeric(val))
stop( i %>% paste('must be Numeric') )
if(val <= 0)
stop( i %>% paste('must be positive') )
}
if(min_density > 1)
stop('min_density must be between 0 and 1')
func = function(las){
if(min(las$Z) < 0)
message("points with Z below 0 will be ignored")
if(min(las$Z) > 5)
message("point cloud doesn't look normalized (Z values too high) - check ?tlsNormalize")
survey_points = if(las %>% hasField('Classification')){
las$Classification != 2
}else{
rep(T, las@data %>% nrow)
}
if(!hasField(las, 'TreeID')){
message('no TreeID field found with tree_points signature: performing single stem point classification')
results = houghStemPoints(las2xyz(las)[survey_points,], h_base[1], h_base[2], h_step, max_d/2, pixel_size, min_density, min_votes)
}else{
message('performing point classification on multiple stems')
survey_points = survey_points & las$TreeID > 0
results = houghStemPlot(las2xyz(las)[survey_points,], las@data$TreeID[survey_points], h_base[1], h_base[2], h_step, max_d/2, pixel_size, min_density, min_votes)
}
las@data$Stem = F
las@data$Stem[survey_points] = results$Stem
las@data$Segment = 0
las@data$Segment[survey_points] = results$Segment
las@data$Radius = 0
las@data$Radius[survey_points] = results$Radius
las@data$Votes = 0
las@data$Votes[survey_points] = results$Votes
las = cleanFields(las, c('Radius', 'Votes'))
return(las)
}
func %<>% setAttribute('stem_pts_mtd')
return(func)
}
#' Stem denoising algorithm: KNN eigen decomposition + point normals intersections voting
#' @description This function is meant to be used inside \code{\link{stemPoints}}. It filters points based on their nearest neighborhood geometries (check \code{\link{fastPointMetrics}}) and assign them to stem patches if reaching a voxel with enough votes.
#' @template param-h_step
#' @template param-max-curvature
#' @template param-max-verticality
#' @param voxel_spacing \code{numeric} - voxel (or pixel) spacing for counting point normals intersections.
#' @template param-max-d
#' @template param-votes-weight
#' @template param-v3d
#' @template section-eigen-decomposition
#' @template section-normals-voting
#' @template reference-liang
#' @template reference-olofsson-2016
#' @export
stm.eigen.knn = function(h_step = .5, max_curvature = .1, max_verticality = 10, voxel_spacing = .025, max_d = .5, votes_weight = .2, v3d = FALSE){
params = list(
h_step = h_step,
max_curvature = max_curvature,
max_verticality = max_verticality,
voxel_spacing = voxel_spacing,
max_d = max_d,
votes_weight = votes_weight
)
for(i in names(params)){
val = params[[i]]
if(length(val) != 1)
stop( i %>% paste('must be of length 1') )
if(!is.numeric(val))
stop( i %>% paste('must be Numeric') )
if(val <= 0)
stop( i %>% paste('must be positive') )
}
if(max_curvature > 1) stop('max_curvature must be a number between 0 and 1')
if(max_verticality > 180) stop('max_verticality must be a number between 0 and 180')
if(votes_weight > 1) stop('votes_weight must be a number between 0 and 1')
func = function(las){
las@data$PointID = 1:nrow(las@data)
mtrlst = c('N', 'Curvature', 'Verticality', 'EigenVector13', 'EigenVector23', 'EigenVector33')
# if(hasField(las, 'TreeID')){
# las@data = las@data %>% split(las@data$TreeID) %>% lapply(LAS) %>%
# lapply(fastPointMetrics, method = ptm.knn(), metrics_list=mtrlst) %>%
# lapply(function(x) x@data) %>% do.call(what=rbind) %>% as.data.table
# }
check_point_metrics = mtrlst %>% sapply(function(x) hasField(las, x)) %>% as.logical %>% all
if(!check_point_metrics){
message('Calculating knn fastPointMetrics')
las = fastPointMetrics(las, ptm.knn(), mtrlst)
}
las@data$Stem = with(las@data, N > 3 & Curvature < max_curvature & abs(Verticality - 90) < max_verticality)
if(hasField(las, 'Classification')){
las@data$Stem = las@data$Stem & las@data$Classification != 2
}
stem_seg = seq(0, max(las$Z)+h_step, h_step)
las@data$Segment = cut(las$Z, stem_seg, include_lowest=T, right=F, ordered_result=T) %>% as.integer
las@data$Segment[las@data$Z < 0] = 0
if(hasField(las, 'TreeID')){
las@data$Stem = las@data$Stem & las@data$TreeID > 0
points = las@data[Stem & order(TreeID, Segment, PointID), .(TreeID, Segment, PointID, X, Y, Z, EigenVector13, EigenVector23, EigenVector33)]
tds = points$TreeID
sgs = points$Segment
ids = points$PointID
votes = points[,-c(1:3)] %>% as.matrix %>% plotEigenHough(ids, tds, sgs, voxel_spacing, max_d/2, !v3d, F) %>%
lapply(function(x) x %>% do.call(what=cbind)) %>% do.call(what=rbind) %>% as.data.table
colnames(votes) = c('Votes','Radius','PointID', 'Segment', 'TreeID')
}else{
points = las@data[Stem & order(Segment, PointID), .(Segment, PointID, X, Y, Z, EigenVector13, EigenVector23, EigenVector33)]
sgs = points$Segment
ids = points$PointID
votes = points[,-c(1:2)] %>% as.matrix %>% treeEigenHough(ids, sgs, voxel_spacing, max_d/2, !v3d, F) %>%
lapply(function(x) x %>% do.call(what=cbind)) %>% do.call(what=rbind) %>% as.data.table
colnames(votes) = c('Votes','Radius','PointID', 'Segment')
}
keepcols = !(colnames(las@data) %in% c('Votes', 'Radius', 'MaxVotes'))
keepcols = colnames(las@data)[keepcols]
las@data = las@data[,keepcols,with=F]
las@data = merge(las@data, votes[,.(PointID, Votes, Radius)], by='PointID', sort=F, all.x=T)
las@data[!las@data$Stem, c('Votes', 'Radius')] = 0
by_cols = 'Segment'
if(hasField(las, 'TreeID')){
max_votes = las@data[,.(MaxVotes = max(Votes)), by=TreeID]
las@data = merge(las@data, max_votes, by='TreeID', sort=FALSE, all.x=TRUE)
las@data$VotesWeight = las@data$Votes / las@data$MaxVotes
las@data$MaxVotes = NULL
by_cols %<>% c('TreeID')
}else{
las@data$VotesWeight = las@data$Votes / max(las@data$Votes)
}
las@data$Stem = las@data$VotesWeight > votes_weight
radii = las@data[Radius > 0 & Stem, .(med_rad = median(Radius)), by=by_cols]
las@data = merge(las@data, radii, by=by_cols, sort=FALSE, all.x=TRUE)
las@data$Radius = las@data$med_rad
las@data$med_rad = NULL
las@data[Stem == F]$Radius = 0
las = cleanFields(las, c('Votes', 'VotesWeight', 'Radius'))
return(las)
}
func %<>% setAttribute('stem_pts_mtd')
return(func)
}
#' Stem denoising algorithm: Voxel eigen decomposition + point normals intersections voting
#' @description This function is meant to be used inside \code{\link{stemPoints}}. It filters points based on their voxel geometries (check \code{\link{fastPointMetrics}}) and assign them to stem patches if reaching a voxel with enough votes.
#' @template param-h_step
#' @template param-max-curvature
#' @template param-max-verticality
#' @template param-voxel-spacing
#' @template param-max-d
#' @template param-votes-weight
#' @template param-v3d
#' @template section-eigen-decomposition
#' @template section-normals-voting
#' @template reference-liang
#' @template reference-olofsson-2016
#' @export
stm.eigen.voxel = function(h_step = .5, max_curvature = .1, max_verticality = 10, voxel_spacing = .025, max_d = .5, votes_weight = .2, v3d = FALSE){
params = list(
h_step = h_step,
max_curvature = max_curvature,
max_verticality = max_verticality,
voxel_spacing = voxel_spacing,
max_d = max_d,
votes_weight = votes_weight
)
for(i in names(params)){
val = params[[i]]
if(length(val) != 1)
stop( i %>% paste('must be of length 1') )
if(!is.numeric(val))
stop( i %>% paste('must be Numeric') )
if(val <= 0)
stop( i %>% paste('must be positive') )
}
if(max_curvature > 1) stop('max_curvature must be a number between 0 and 1')
if(max_verticality > 180) stop('max_verticality must be a number between 0 and 180')
if(votes_weight > 1) stop('votes_weight must be a number between 0 and 1')
func = function(las){
mtrlst = c('N', 'Curvature', 'Verticality', 'EigenVector13', 'EigenVector23', 'EigenVector33', 'VoxelID')
check_point_metrics = mtrlst %>% sapply(function(x) hasField(las, x)) %>% as.logical %>% all
if(!check_point_metrics){
message('Calculating voxel fastPointMetrics')
las = fastPointMetrics(las, ptm.voxel(voxel_spacing), mtrlst)
}
las@data$Stem = with(las@data, Classification != 2 & N > 3 & Curvature < max_curvature & abs(Verticality - 90) < max_verticality)
stem_seg = seq(0, max(las$Z)+h_step, h_step)
las@data$Segment = cut(las$Z, stem_seg, include_lowest=T, right=F, ordered_result=T) %>% as.integer
las@data$Segment[las@data$Z < 0] = 0
if(hasField(las, 'TreeID')){
voxels = las@data[Stem & order(TreeID, Segment), .(TreeID = mean(TreeID), Segment = mean(Segment), X = mean(X), Y = mean(Y), Z = mean(Z), EigenVector13 = mean(EigenVector13), EigenVector23 = mean(EigenVector23), EigenVector33 = mean(EigenVector33)), by = VoxelID]
tds = voxels$TreeID
sgs = voxels$Segment
ids = voxels$VoxelID
votes = voxels[,-c(1:3)] %>% as.matrix %>% plotEigenHough(ids, tds, sgs, voxel_spacing, max_d/2, !v3d, F) %>%
lapply(function(x) x %>% do.call(what=cbind)) %>% do.call(what=rbind) %>% as.data.table
colnames(votes) = c('Votes','Radius','VoxelID', 'Segment', 'TreeID')
}else{
voxels = las@data[Stem & order(Segment), .(Segment = mean(Segment), X = mean(X), Y = mean(Y), Z = mean(Z), EigenVector13 = mean(EigenVector13), EigenVector23 = mean(EigenVector23), EigenVector33 = mean(EigenVector33)), by = VoxelID]
sgs = voxels$Segment
ids = voxels$VoxelID
votes = voxels[,-c(1:2)] %>% as.matrix %>% treeEigenHough(ids, sgs, voxel_spacing, max_d/2, !v3d, F) %>%
lapply(function(x) x %>% do.call(what=cbind)) %>% do.call(what=rbind) %>% as.data.table
colnames(votes) = c('Votes','Radius','VoxelID', 'Segment')
}
keepcols = !(colnames(las@data) %in% c('Votes', 'Radius', 'MaxVotes'))
keepcols = colnames(las@data)[keepcols]
las@data = las@data[,keepcols,with=F]
las@data = merge(las@data, votes[,.(VoxelID, Votes, Radius)], by='VoxelID', sort=F, all.x=T)
las@data[!las@data$Stem, c('Votes', 'Radius')] = 0
by_cols = 'Segment'
if(hasField(las, 'TreeID')){
max_votes = las@data[,.(MaxVotes = max(Votes)), by=TreeID]
las@data = merge(las@data, max_votes, by='TreeID', sort=F, all.x=T)
las@data$VotesWeight = las@data$Votes / las@data$MaxVotes
las@data$MaxVotes = NULL
by_cols %<>% c('TreeID')
}else{
las@data$VotesWeight = las@data$Votes / max(las@data$Votes)
}
las@data$Stem = las@data$VotesWeight > votes_weight
radii = las@data[Radius > 0 & Stem, .(med_rad = median(Radius)), by=by_cols]
las@data = merge(las@data, radii, by=by_cols, sort=FALSE, all.x=TRUE)
las@data$Radius = las@data$med_rad
las@data$med_rad = NULL
las@data[Stem == F]$Radius = 0
las = cleanFields(las, c('Votes', 'VotesWeight', 'Radius'))
return(las)
}
func %<>% setAttribute('stem_pts_mtd')
return(func)
}
Try the TreeLS package in your browser
Any scripts or data that you put into this service are public.
TreeLS documentation built on Aug. 26, 2020, 5:14 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions stm.eigen.knn stm.hough
Documented in stm.eigen.knn stm.hough
# ===============================================================================
#
# Developers:
#
# Tiago de Conto - tdc.florestal@gmail.com - https://github.com/tiagodc/
#
# COPYRIGHT: Tiago de Conto, 2020
#
# This piece of software is open and free to use, redistribution and modifications
# should be done in accordance to the GNU General Public License >= 3
# Use this software as you wish, but no warranty is provided whatsoever. For any
# comments or questions on TreeLS, please contact the developer (prefereably through my github account)
#
# If publishing any work/study/research that used the tools in TreeLS,
# please don't forget to cite the proper sources!
#
# Enjoy!
#
# ===============================================================================
#' Stem denoising algorithm: Hough Transform
#' @description This function is meant to be used inside \code{\link{stemPoints}}. It applies an adapted version of the Hough Transform for circle search. Mode details are given in the sections below.
#' @template param-h_step
#' @template param-max-d
#' @template param-hbase
#' @template param-pixel-size
#' @template param-min-density
#' @template param-min-votes
#' @section \code{LAS@data} Special Fields:
#'
#' Meaninful new fields in the output:
#'
#' \itemize{
#' \item \code{Stem}: \code{TRUE} for stem points
#' \item \code{Segment}: stem segment number (from bottom to top and nested with TreeID)
#' \item \code{Radius}: approximate radius of the point's stem segment estimated by the Hough Transform - always a multiple of the \code{pixel_size}
#' \item \code{Votes}: votes received by the stem segment's center through the Hough Transform
#' }
#'
#' @template section-hough-transform
#' @template reference-olofsson
#' @template reference-thesis
#' @export
stm.hough = function(h_step=0.5, max_d=0.5, h_base = c(1,2.5), pixel_size=0.025, min_density=0.1, min_votes=3){
if(length(h_base) != 2)
stop('h_base must be a numeric vector of length 2')
if(diff(h_base) <= 0)
stop('h_base[2] must be larger than h_base[1]')
params = list(
h_step = h_step,
max_d = max_d,
pixel_size = pixel_size,
min_density = min_density,
min_votes = min_votes
)
for(i in names(params)){
val = params[[i]]
if(length(val) != 1)
stop( i %>% paste('must be of length 1') )
if(!is.numeric(val))
stop( i %>% paste('must be Numeric') )
if(val <= 0)
stop( i %>% paste('must be positive') )
}
if(min_density > 1)
stop('min_density must be between 0 and 1')
func = function(las){
if(min(las$Z) < 0)
message("points with Z below 0 will be ignored")
if(min(las$Z) > 5)
message("point cloud doesn't look normalized (Z values too high) - check ?tlsNormalize")
survey_points = if(las %>% hasField('Classification')){
las$Classification != 2
}else{
rep(T, las@data %>% nrow)
}
if(!hasField(las, 'TreeID')){
message('no TreeID field found with tree_points signature: performing single stem point classification')
results = houghStemPoints(las2xyz(las)[survey_points,], h_base[1], h_base[2], h_step, max_d/2, pixel_size, min_density, min_votes)
}else{
message('performing point classification on multiple stems')
survey_points = survey_points & las$TreeID > 0
results = houghStemPlot(las2xyz(las)[survey_points,], las@data$TreeID[survey_points], h_base[1], h_base[2], h_step, max_d/2, pixel_size, min_density, min_votes)
}
las@data$Stem = F
las@data$Stem[survey_points] = results$Stem
las@data$Segment = 0
las@data$Segment[survey_points] = results$Segment
las@data$Radius = 0
las@data$Radius[survey_points] = results$Radius
las@data$Votes = 0
las@data$Votes[survey_points] = results$Votes
las = cleanFields(las, c('Radius', 'Votes'))
return(las)
}
func %<>% setAttribute('stem_pts_mtd')
return(func)
}
#' Stem denoising algorithm: KNN eigen decomposition + point normals intersections voting
#' @description This function is meant to be used inside \code{\link{stemPoints}}. It filters points based on their nearest neighborhood geometries (check \code{\link{fastPointMetrics}}) and assign them to stem patches if reaching a voxel with enough votes.
#' @template param-h_step
#' @template param-max-curvature
#' @template param-max-verticality
#' @param voxel_spacing \code{numeric} - voxel (or pixel) spacing for counting point normals intersections.
#' @template param-max-d
#' @template param-votes-weight
#' @template param-v3d
#' @template section-eigen-decomposition
#' @template section-normals-voting
#' @template reference-liang
#' @template reference-olofsson-2016
#' @export
stm.eigen.knn = function(h_step = .5, max_curvature = .1, max_verticality = 10, voxel_spacing = .025, max_d = .5, votes_weight = .2, v3d = FALSE){
params = list(
h_step = h_step,
max_curvature = max_curvature,
max_verticality = max_verticality,
voxel_spacing = voxel_spacing,
max_d = max_d,
votes_weight = votes_weight
)
for(i in names(params)){
val = params[[i]]
if(length(val) != 1)
stop( i %>% paste('must be of length 1') )
if(!is.numeric(val))
stop( i %>% paste('must be Numeric') )
if(val <= 0)
stop( i %>% paste('must be positive') )
}
if(max_curvature > 1) stop('max_curvature must be a number between 0 and 1')
if(max_verticality > 180) stop('max_verticality must be a number between 0 and 180')
if(votes_weight > 1) stop('votes_weight must be a number between 0 and 1')
func = function(las){
las@data$PointID = 1:nrow(las@data)
mtrlst = c('N', 'Curvature', 'Verticality', 'EigenVector13', 'EigenVector23', 'EigenVector33')
# if(hasField(las, 'TreeID')){
# las@data = las@data %>% split(las@data$TreeID) %>% lapply(LAS) %>%
# lapply(fastPointMetrics, method = ptm.knn(), metrics_list=mtrlst) %>%
# lapply(function(x) x@data) %>% do.call(what=rbind) %>% as.data.table
# }
check_point_metrics = mtrlst %>% sapply(function(x) hasField(las, x)) %>% as.logical %>% all
if(!check_point_metrics){
message('Calculating knn fastPointMetrics')
las = fastPointMetrics(las, ptm.knn(), mtrlst)
}
las@data$Stem = with(las@data, N > 3 & Curvature < max_curvature & abs(Verticality - 90) < max_verticality)
if(hasField(las, 'Classification')){
las@data$Stem = las@data$Stem & las@data$Classification != 2
}
stem_seg = seq(0, max(las$Z)+h_step, h_step)
las@data$Segment = cut(las$Z, stem_seg, include_lowest=T, right=F, ordered_result=T) %>% as.integer
las@data$Segment[las@data$Z < 0] = 0
if(hasField(las, 'TreeID')){
las@data$Stem = las@data$Stem & las@data$TreeID > 0
points = las@data[Stem & order(TreeID, Segment, PointID), .(TreeID, Segment, PointID, X, Y, Z, EigenVector13, EigenVector23, EigenVector33)]
tds = points$TreeID
sgs = points$Segment
ids = points$PointID
votes = points[,-c(1:3)] %>% as.matrix %>% plotEigenHough(ids, tds, sgs, voxel_spacing, max_d/2, !v3d, F) %>%
lapply(function(x) x %>% do.call(what=cbind)) %>% do.call(what=rbind) %>% as.data.table
colnames(votes) = c('Votes','Radius','PointID', 'Segment', 'TreeID')
}else{
points = las@data[Stem & order(Segment, PointID), .(Segment, PointID, X, Y, Z, EigenVector13, EigenVector23, EigenVector33)]
sgs = points$Segment
ids = points$PointID
votes = points[,-c(1:2)] %>% as.matrix %>% treeEigenHough(ids, sgs, voxel_spacing, max_d/2, !v3d, F) %>%
lapply(function(x) x %>% do.call(what=cbind)) %>% do.call(what=rbind) %>% as.data.table
colnames(votes) = c('Votes','Radius','PointID', 'Segment')
}
keepcols = !(colnames(las@data) %in% c('Votes', 'Radius', 'MaxVotes'))
keepcols = colnames(las@data)[keepcols]
las@data = las@data[,keepcols,with=F]
las@data = merge(las@data, votes[,.(PointID, Votes, Radius)], by='PointID', sort=F, all.x=T)
las@data[!las@data$Stem, c('Votes', 'Radius')] = 0
by_cols = 'Segment'
if(hasField(las, 'TreeID')){
max_votes = las@data[,.(MaxVotes = max(Votes)), by=TreeID]
las@data = merge(las@data, max_votes, by='TreeID', sort=FALSE, all.x=TRUE)
las@data$VotesWeight = las@data$Votes / las@data$MaxVotes
las@data$MaxVotes = NULL
by_cols %<>% c('TreeID')
}else{
las@data$VotesWeight = las@data$Votes / max(las@data$Votes)
}
las@data$Stem = las@data$VotesWeight > votes_weight
radii = las@data[Radius > 0 & Stem, .(med_rad = median(Radius)), by=by_cols]
las@data = merge(las@data, radii, by=by_cols, sort=FALSE, all.x=TRUE)
las@data$Radius = las@data$med_rad
las@data$med_rad = NULL
las@data[Stem == F]$Radius = 0
las = cleanFields(las, c('Votes', 'VotesWeight', 'Radius'))
return(las)
}
func %<>% setAttribute('stem_pts_mtd')
return(func)
}
#' Stem denoising algorithm: Voxel eigen decomposition + point normals intersections voting
#' @description This function is meant to be used inside \code{\link{stemPoints}}. It filters points based on their voxel geometries (check \code{\link{fastPointMetrics}}) and assign them to stem patches if reaching a voxel with enough votes.
#' @template param-h_step
#' @template param-max-curvature
#' @template param-max-verticality
#' @template param-voxel-spacing
#' @template param-max-d
#' @template param-votes-weight
#' @template param-v3d
#' @template section-eigen-decomposition
#' @template section-normals-voting
#' @template reference-liang
#' @template reference-olofsson-2016
#' @export
stm.eigen.voxel = function(h_step = .5, max_curvature = .1, max_verticality = 10, voxel_spacing = .025, max_d = .5, votes_weight = .2, v3d = FALSE){
params = list(
h_step = h_step,
max_curvature = max_curvature,
max_verticality = max_verticality,
voxel_spacing = voxel_spacing,
max_d = max_d,
votes_weight = votes_weight
)
for(i in names(params)){
val = params[[i]]
if(length(val) != 1)
stop( i %>% paste('must be of length 1') )
if(!is.numeric(val))
stop( i %>% paste('must be Numeric') )
if(val <= 0)
stop( i %>% paste('must be positive') )
}
if(max_curvature > 1) stop('max_curvature must be a number between 0 and 1')
if(max_verticality > 180) stop('max_verticality must be a number between 0 and 180')
if(votes_weight > 1) stop('votes_weight must be a number between 0 and 1')
func = function(las){
mtrlst = c('N', 'Curvature', 'Verticality', 'EigenVector13', 'EigenVector23', 'EigenVector33', 'VoxelID')
check_point_metrics = mtrlst %>% sapply(function(x) hasField(las, x)) %>% as.logical %>% all
if(!check_point_metrics){
message('Calculating voxel fastPointMetrics')
las = fastPointMetrics(las, ptm.voxel(voxel_spacing), mtrlst)
}
las@data$Stem = with(las@data, Classification != 2 & N > 3 & Curvature < max_curvature & abs(Verticality - 90) < max_verticality)
stem_seg = seq(0, max(las$Z)+h_step, h_step)
las@data$Segment = cut(las$Z, stem_seg, include_lowest=T, right=F, ordered_result=T) %>% as.integer
las@data$Segment[las@data$Z < 0] = 0
if(hasField(las, 'TreeID')){
voxels = las@data[Stem & order(TreeID, Segment), .(TreeID = mean(TreeID), Segment = mean(Segment), X = mean(X), Y = mean(Y), Z = mean(Z), EigenVector13 = mean(EigenVector13), EigenVector23 = mean(EigenVector23), EigenVector33 = mean(EigenVector33)), by = VoxelID]
tds = voxels$TreeID
sgs = voxels$Segment
ids = voxels$VoxelID
votes = voxels[,-c(1:3)] %>% as.matrix %>% plotEigenHough(ids, tds, sgs, voxel_spacing, max_d/2, !v3d, F) %>%
lapply(function(x) x %>% do.call(what=cbind)) %>% do.call(what=rbind) %>% as.data.table
colnames(votes) = c('Votes','Radius','VoxelID', 'Segment', 'TreeID')
}else{
voxels = las@data[Stem & order(Segment), .(Segment = mean(Segment), X = mean(X), Y = mean(Y), Z = mean(Z), EigenVector13 = mean(EigenVector13), EigenVector23 = mean(EigenVector23), EigenVector33 = mean(EigenVector33)), by = VoxelID]
sgs = voxels$Segment
ids = voxels$VoxelID
votes = voxels[,-c(1:2)] %>% as.matrix %>% treeEigenHough(ids, sgs, voxel_spacing, max_d/2, !v3d, F) %>%
lapply(function(x) x %>% do.call(what=cbind)) %>% do.call(what=rbind) %>% as.data.table
colnames(votes) = c('Votes','Radius','VoxelID', 'Segment')
}
keepcols = !(colnames(las@data) %in% c('Votes', 'Radius', 'MaxVotes'))
keepcols = colnames(las@data)[keepcols]
las@data = las@data[,keepcols,with=F]
las@data = merge(las@data, votes[,.(VoxelID, Votes, Radius)], by='VoxelID', sort=F, all.x=T)
las@data[!las@data$Stem, c('Votes', 'Radius')] = 0
by_cols = 'Segment'
if(hasField(las, 'TreeID')){
max_votes = las@data[,.(MaxVotes = max(Votes)), by=TreeID]
las@data = merge(las@data, max_votes, by='TreeID', sort=F, all.x=T)
las@data$VotesWeight = las@data$Votes / las@data$MaxVotes
las@data$MaxVotes = NULL
by_cols %<>% c('TreeID')
}else{
las@data$VotesWeight = las@data$Votes / max(las@data$Votes)
}
las@data$Stem = las@data$VotesWeight > votes_weight
radii = las@data[Radius > 0 & Stem, .(med_rad = median(Radius)), by=by_cols]
las@data = merge(las@data, radii, by=by_cols, sort=FALSE, all.x=TRUE)
las@data$Radius = las@data$med_rad
las@data$med_rad = NULL
las@data[Stem == F]$Radius = 0
las = cleanFields(las, c('Votes', 'VotesWeight', 'Radius'))
return(las)
}
func %<>% setAttribute('stem_pts_mtd')
return(func)
}
Try the TreeLS package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.