In tommywhitney/mobile_lidar: What the Package Does (One Line, Title Case)
library(tidyverse)
library(scatterplot3d)
library(rgl)
library(rospca)
library(geometry)
library(spanner)
library(dbscan)
library(ggforce)
library(cowplot)
library(Rcpp)
library(mobileLidar)
library(microbenchmark)
# las <- readTLSLAS(file.choose())
# las <- readTLSLAS("C:/Users/matt/Desktop/DensePatchA.laz")
# Rudd tank burn only north
# las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_BO_N_zebcam_ALSregistered.laz", filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431804.383088774 3899195.49902686 30 (x y radius)")
# rudd tank burn only south east
# las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_BO_SE_zebcam_ALSregistered.laz", filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431838.287584286 3899147.61210088 30 (x y radius)")
# # rudd tank burn only south west
las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_BO_SW_zebcam_ALSregistered.laz", select = "xyz", filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431790.667864794 3899145.91081599 30 (x y radius)")
# control north
# las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_CO_N_zebcam_ALSregistered.laz",
# select = "xyz",
# filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431554.241897831 3899192.99624304 30 (x y radius)")
# # control south east
# las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_CO_SE_zebcam_ALSregistered.laz",
# select = "xyz",
# filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431595.891142449 3899156.16591304 30 (x y radius)")
# control south west
# las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_CO_SW_zebcam_ALSregistered.laz",
# select = "xyz",
# filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431545.984681478 3899149.36050936 30 (x y radius)")
# treat and burn north east
# las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_TB_NE_zebcam_ALSregistered.laz",
# select = "xyz",
# filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431940.528400044 3899570.10647836 30 (x y radius)")
# treat and burn north west
# las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_TB_NW_zebcam_ALSregistered.laz",
# select = "xyz",
# filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431883.255561657 3899552.23327632 30 (x y radius)")
# # treat and burn south
# las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_TB_S_zebcam_ALSregistered.laz",
# select = "xyz",
# filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431922.524931294 3899521.26763958 30 (x y radius)")
# # treat north
# las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_TO_N_zebcam_ALSregistered.laz",
# select = "xyz",
# filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431549.858895073 3899580.94101544 30 (x y radius)")
# # treat south east
# las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_TO_SE_zebcam_ALSregistered.laz",
# select = "xyz",
# filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431572.236974267 3899519.5420316 30 (x y radius)")
# # treat south west
# las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_TO_SW_zebcam_ALSregistered.laz",
# select = "xyz",
# filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431506.243392164 3899537.98026232 30 (x y radius)")
las <- classify_ground(las, csf(sloop_smooth = FALSE,
class_threshold = 0.5,
cloth_resolution = 0.5,
rigidness = 1L,
iterations = 500L,
time_step = 0.65))
las <- normalize_height(las, tin())
las <- classify_noise(las, ivf(6,1))
las <- filter_poi(las, Classification != LASNOISE)
plot(las, axis = T)
las_slice <- filter_poi(las, Z>=0.62, Z<=2.12)
eigens <- spanner::eigen_metrics(las_slice, radius = 0.33, ncpu = 8)
pt_den <- spanner:::C_count_in_sphere(las_slice, radius = 0.33, ncpu = 8)
las_slice@data<-cbind(las_slice@data, eigens)
las_slice@data<-cbind(las_slice@data, pt_den)
# adjust this Z slice if the trees width is close to the height of the slice
las_slice <- filter_poi(las_slice, Z>=0.62 ,Z<= 2.12)
clust <- dbscan::dbscan(las_slice@data[,c("X","Y","Z", "eSum","Verticality")], eps = 0.25, minPts = 100)
las_slice@data$treeID<-clust$cluster
plot(las_slice, color="treeID", axis = T)
fit_df <- las_slice_circle_fitting( las_slice, 80, 0.05, 0.6)
fit_df <- fit_df[fit_df$Radius <= 2, ]
offsets <- plot(las_slice, color="treeID", axis = T)
spheres3d( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = 1.37, r = fit_df[,3], alpha = .7)
rgl.texts( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = fit_df[,3]+1.37 +.1, text = paste( "x: ",fit_df[,1]))
rgl.texts( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = fit_df[,3]+1.37 +.2, text = paste( "y: ", fit_df[,2] ))
rgl.texts( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = fit_df[,3]+1.37 +.3, text = paste( "diameter: ", fit_df[,3]*2))
rgl.texts( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = fit_df[,3]+1.37 +.4, text = paste( "mean squared error: ", fit_df[,4] ))
rgl.texts( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = fit_df[,3]+1.37 +.5, text = paste( "inclusion: ", fit_df[,5] ))
rgl.texts( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = fit_df[,3]+1.37 +.6, text = paste( "tree id: ", fit_df[,6] ))
rgl.texts( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = fit_df[,3]+1.37 +.7, text = paste( "lean: ", fit_df[,7] ))
test_tree = 1
tree_slice <- filter_poi(las_slice, treeID == test_tree)
# call las_slice_circle_fitting on tree slice to get a new fit with given parameters
new_fit <- las_slice_circle_fitting( tree_slice, 80, 0.01, 0.6 )
fit_df <- rbind(fit_df, new_fit)
plot(tree_slice, color="treeID", axis = T)
# assign all points to a matrix
M <- cbind( tree_slice$X, tree_slice$Y, tree_slice$Z, tree_slice$treeID)
# select a tree from the matrix by id
# tree_pts <- M[ M[ ,4] == i, ]
# assign tmp matrix the points for current tree id
tmp <- M[ , 1:3]
# preform robust pca
resR <- robpca(tmp, k=3, ndir = 5000)
# filter tmp to around dbh
# tmp <- tmp[ tmp[ ,3] >= 1.27, ]
# tmp <- tmp[ tmp[ ,3] <= 1.47, ]
# use robust pca to reduce dimensions and fit points to a plane orthogonal to the tree lean
# get mean and principal components of points
# note here pc1 will be the greatest value of the principal componets
# with pc1 being the trees vertical componet, and pc2 and pc3 being the horizontal components
# if the tree has a greater width than the height of the slice taken, what we assume to be pc1 as
# the height will be one of the horozontal componets
center <- resR$center
pc1 <- resR$loadings[ , 1]
pc2 <- resR$loadings[ , 2]
pc3 <- resR$loadings[ , 3]
# bind the lower two principal components into a matrix
pc <- t(cbind( pc3, pc2))
# subtract means from each point
tmp_sub_means <- t( sweep( tmp, 2, center ))
# init empty matrix
pts_proj <- matrix( 0,nrow = 2)
# loop through all points passing them through matrix
for( i in 1:ncol(tmp_sub_means))
{
a <- pc %*% tmp_sub_means[ , i ]
pts_proj <- cbind(pts_proj, a)
}
# for some reason when removing these na values the loop does not end
# remove first row of matrix, contains na values
# pts_proj <-pts_proj[ -1,]
# transpose matrix to long format
pts_proj <- t(pts_proj)
pts_proj <- pts_proj[ -c(1),]
df <- as.data.frame( pts_proj)
best_fit <- ransac_circle_fit( pts_proj, 10000, 0.1, 0.8)
format(best_fit, scientific = F)
center
xyz_center <- center + best_fit[1]%*%pc3 + best_fit[2]%*%pc2
plot(tree_slice, color="treeID")
# spheres3d(x = xyz_center[1], y = xyz_center[2], z = 1.37, radius = best_fit[3])
ggplot( df, aes( x = pc3, y = pc2))+
geom_point()+
ggtitle("tree points projected")+
geom_circle( aes( x0 = best_fit[1], y0 = best_fit[2], r = best_fit[3], color = "red"), inherit.aes = F )+
coord_fixed()+
theme( legend.position = "none")+
geom_point(aes(x=best_fit[1], y=best_fit[2]), colour="blue")+
coord_fixed()
offsets <- plot(tree_slice, color="treeID")
spheres3d(x = xyz_center[,1]-offsets[1], y = xyz_center[,2]-offsets[2], z = 1.37, r = fit_df[,3], alpha = .7)
# rgl.texts( x = xyz_center[,1]-offsets[1], y = xyz_center[,2]-offsets[2], z = 2.3, text = paste( "diameter: ", fit_df[,3]*2))
# rgl.texts( x = xyz_center[,1]-offsets[1], y = xyz_center[,2]-offsets[2], z = 2.15, text = paste( "x: ",fit_df[,1]))
# rgl.texts( x = xyz_center[,1]-offsets[1], y = xyz_center[,2]-offsets[2], z = 2, text = paste( "y: ", fit_df[,2] ))
open3d()
theta <- seq(0, 2*pi, len = 25)
knot <- cylinder3d(
center = cbind(
sin(theta) + 2*sin(2*theta),
2*sin(3*theta),
cos(theta) - 2*cos(2*theta)),
e1 = cbind(
cos(theta) + 4*cos(2*theta),
6*cos(3*theta),
sin(theta) + 4*sin(2*theta)),
radius = 0.8,
closed = TRUE)
# assign all points to a matrix
M <- cbind( tree_slice$X, tree_slice$Y, tree_slice$Z, tree_slice$treeID)
# assign tmp matrix the points for current tree id
tmp <- M[ , 1:3]
# preform robust pca
resR <- robpca(tmp, k=3, ndir = 5000)
# filter tmp to around dbh
# tmp <- tmp[ tmp[ ,3] >= 1.27, ]
# tmp <- tmp[ tmp[ ,3] <= 1.47, ]
# use robust pca to reduce dimensions and fit points to a plane orthogonal to the tree lean
# get mean and principal components of points
# note here pc1 will be the greatest value of the principal componets
# with pc1 being the trees vertical componet, and pc2 and pc3 being the horizontal components
# if the tree has a greater width than the height of the slice taken, what we assume to be pc1 as
# the height will be one of the horozontal componets
center <- resR$center
pc1 <- resR$loadings[ , 1]
pc2 <- resR$loadings[ , 2]
pc3 <- resR$loadings[ , 3]
# bind the lower two principal components into a matrix
pc <- t(cbind( pc3, pc2))
# subtract means from each point
tmp_sub_means <- t( sweep( tmp, 2, center ))
# init empty matrix
pts_proj <- matrix( 0,nrow = 2)
# loop through all points passing them through matrix
for( i in 1:ncol(tmp_sub_means))
{
a <- pc %*% tmp_sub_means[ , i ]
pts_proj <- cbind(pts_proj, a)
}
# for some reason when removing these na values the loop does not end
# remove first row of matrix, contains na values
# pts_proj <-pts_proj[ -1,]
# transpose matrix to long format
pts_proj <- t(pts_proj)
pts_proj <- pts_proj[ -c(1),]
df <- as.data.frame( pts_proj)
ggplot( df, aes( x = pc3, y = pc2))+
geom_point()+
ggtitle("tree points projected")+
geom_circle( aes( x0 = fit_df[test_tree,1], y0 = fit_df[test_tree,2], r = fit_df[test_tree,3], color = "red"), inherit.aes = F )+
coord_fixed()+
theme( legend.position = "none")+
coord_fixed()
fit_df[test_tree,]
# assign all points to a matrix
M <- cbind( las_slice$X, las_slice$Y, las_slice$Z, las_slice$treeID)
# select a tree from the matrix by id
tree_pts <- M[ M[ ,4] == i, ]
# assign tmp matrix the points for current tree id
tmp <- tree_pts[ , 1:3]
# preform robust pca
resR <- robpca(tmp, k=3, ndir = 5000)
# filter tmp to around dbh
# tmp <- tmp[ tmp[ ,3] >= 1.27, ]
# tmp <- tmp[ tmp[ ,3] <= 1.47, ]
# use robust pca to reduce dimensions and fit points to a plane orthogonal to the tree lean
# get mean and principal components of points
# note here pc1 will be the greatest value of the principal componets
# with pc1 being the trees vertical componet, and pc2 and pc3 being the horizontal components
# if the tree has a greater width than the height of the slice taken, what we assume to be pc1 as
# the height will be one of the horozontal componets
center <- resR$center
pc1 <- resR$loadings[ , 1]
pc2 <- resR$loadings[ , 2]
pc3 <- resR$loadings[ , 3]
# bind the lower two principal components into a matrix
pc <- t(cbind( pc3, pc2))
# subtract means from each point
tmp_sub_means <- t( sweep( tmp, 2, center ))
# init empty matrix
pts_proj <- matrix( 0,nrow = 2)
# loop through all points passing them through matrix
for( i in 1:ncol(tmp_sub_means))
{
a <- pc %*% tmp_sub_means[ , i ]
pts_proj <- cbind(pts_proj, a)
}
# for some reason when removing these na values the loop does not end
# remove first row of matrix, contains na values
# pts_proj <-pts_proj[ -1,]
# transpose matrix to long format
pts_proj <- t(pts_proj)
pts_proj <- pts_proj[ -c(1),]
df <- as.data.frame( pts_proj)
best_fit <- ransac_circle_fit( pts_proj, 1000, 0.01, 0.8)
format(best_fit, scientific = F)
best_fit[3]*2
geom_point()+
ggtitle("tree points projected")+
geom_circle( aes( x0 = best_fit[1], y0 = best_fit[2], r = best_fit[3], color = "red"), inherit.aes = F )+
coord_fixed()+
theme( legend.position = "none")+
geom_point(aes(x=best_fit[1], y=best_fit[2]), colour="blue")+
coord_fixed()
rgl.open()
rgl.points(tree_pts[ , 1 ], tree_pts[ , 2], tree_pts[ , 3], color ="white", axis = T)
# select a tree from the matrix by id
tree_pts <- M[ M[ ,4] == 1, ]
# assign tmp matrix the points for current tree id
tmp <- tree_pts[ , 1:3]
# preform robust pca
resR <- robpca(tmp, k=3, ndir = 5000)
# use robust pca to reduce dimensions and fit points to a plane orthogonal to the tree lean
# get mean and principal components of points
# note here pc1 will be the greatest value of the principal componets
# with pc1 being the trees vertical componet, and pc2 and pc3 being the horizontal components
# if the tree has a greater width than the height of the slice taken, what we assume to be pc1 as
# the height will be one of the horozontal componets
center <- resR$center
pc1 <- resR$loadings[ , 1]
pc2 <- resR$loadings[ , 2]
pc3 <- resR$loadings[ , 3]
# bind the lower two principal components into a matrix
pc <- t(cbind( pc3, pc2))(pts_proj)
# subtract means from each point
tmp_sub_means <- t( sweep( tmp, 2, center ))
# init empty matrix
pts_proj <- matrix( nrow = 2)
# loop through all points passing them through matrix
for( i in 1:ncol(tmp_sub_means))
{
a <- pc %*% tmp_sub_means[ , i ]
pts_proj <- cbind(pts_proj, a)
}
# transpose matrix to long format
pts_proj <- t
df <- as.data.frame( pts_proj)
# }
best_fit <- ransac_circle_fit( pts_proj, 10000, 0.001, 0.9)
best_fit
circle fitting
ggplot( df, aes( x = pc3, y = pc2))+
geom_point()+
ggtitle("tree points projected")+
geom_circle( aes( x0 = best_fit[1], y0 = best_fit[2], r = best_fit[3], color = "red"), inherit.aes = F )+
coord_fixed()+
theme( legend.position = "none")
best_fit
best_fit[3]*2
fit at 20cm around dbh
# assign all points to a matrix
M <- cbind( las_slice$X, las_slice$Y, las_slice$Z, las_slice$treeID)
# select a tree from the matrix by id
tree_pts <- M[M[ ,4] == 1,]
# assign tmp matrix the points for current tree id
tmp <- tree_pts[ , 1:3]
# preform robust pca
resR <- robpca(tmp, k=3, ndir = 5000)
# select 10 cm of points around dbh
tmp <- tmp[ tmp[ , 3] >= 1.27 & tmp[ , 3] <= 1.47, ]
# use robust pca to reduce dimensions and fit points to a plane orthogonal to the tree lean
# get mean and principal components of points
# note here pc1 will be the greatest value of the principal componets
# with pc1 being the trees vertical componet, and pc2 and pc3 being the horizontal componets
# if the tree has a greater width than the height of the slice taken what we assume to be pc1 as
# the height will be one of the horozontal componets
center <- resR$center
pc1 <- resR$loadings[ , 1]
pc2 <- resR$loadings[ , 2]
pc3 <- resR$loadings[ , 3]
# bind the lower two principal components into a matrix
pc <- t(cbind( pc3, pc2))
# subtract means from each point
tmp_sub_means <- t( sweep( tmp, 2, center ))
# init empty matrix
pts_proj <- matrix( nrow = 2)
# loop through all points passing them through matrix
for( i in 1:ncol(tmp_sub_means))
{
a <- pc %*% tmp_sub_means[ , i ]
pts_proj <- cbind(pts_proj, a)
}
# transpose matrix to long format
pts_proj <- t(pts_proj)
df <- as.data.frame( pts_proj)
best_fit <- ransac_circle_fit( pts_proj, 10000, 0.001, 0.9)
best_fit
ggplot( df, aes( x = pc2, y = pc3))+
geom_point()+
ggtitle("tree points projected")+
geom_circle( aes( x0 = best_fit[2], y0 = best_fit[1], r = best_fit[3], color = "red"), inherit.aes = F )+
coord_fixed()+
theme( legend.position = "none")
circle fitting
sample <- pts_proj[ sample( nrow(pts_proj), size = 3, replace = F), ]
line_1_y_diff = sample[2, 2] - sample[1, 2]
line_1_x_diff = sample[2, 1] - sample[1, 1]
line_2_y_diff = sample[3, 2] - sample[2, 2]
line_2_x_diff = sample[3, 1] - sample[2, 1]
line_1_slope = line_1_y_diff/line_1_x_diff
line_2_slope = line_2_y_diff/line_2_x_diff
center_x = ( line_1_slope * line_2_slope *
( sample[1, 2] - sample[3, 2]) +
line_2_slope *
( sample[1, 1] + sample[2, 1]) -
line_1_slope *
( sample[2, 1] + sample[3, 1])) /
(2* (line_2_slope-line_1_slope) )
center_y = -1 * (center_x - ( sample[1, 1] + sample[2, 1])/2) /
line_1_slope +
(sample[1, 2]+sample[2, 2])/2
r = sqrt( (center_x - sample[1, 1])^2 + (center_y - sample[1, 2])^2)
center_x
center_y
r
view single iteration circle fit
ggplot( df, aes( x = pc2, y = pc3))+
geom_point()+
ggtitle("tree points projected")+
geom_circle( aes( x0 = center_y, y0 = center_x, r = r, color = "red"), inherit.aes = F )+
coord_fixed()+
theme( legend.position = "none")
# assign all points to a matrix
M <- cbind( las_slice$X, las_slice$Y, las_slice$Z, las_slice$treeID)
# select a tree from the matrix by id
tree_pts <- M[M[ ,4] == 2,]
# assign tmp matrix the points for current tree id
tmp <- tree_pts[ , 1:3]
# preform robust pca
resR <- robpca(tmp, k=3, ndir = 5000)
# select 10 cm of points around dbh
tmp <- tmp[ tmp[ , 3] >= 1.27 & tmp[ , 3] <= 1.47, ]
# use robust pca to reduce dimensions and fit points to a plane orthogonal to the tree lean
# get mean and principal components of points
# note here pc1 will be the greatest value of the principal componets
# with pc1 being the trees vertical componet, and pc2 and pc3 being the horizontal componets
# if the tree has a greater width than the height of the slice taken what we assume to be pc1 as
# the height will be one of the horozontal componets
center <- resR$center
pc1 <- resR$loadings[ , 1]
pc2 <- resR$loadings[ , 2]
pc3 <- resR$loadings[ , 3]
# bind the lower two principal components into a matrix
pc <- t(cbind( pc3, pc2))
# subtract means from each point
tmp_sub_means <- t( sweep( tmp, 2, center ))
# init empty matrix
pts_proj <- matrix( nrow = 2)
# loop through all points passing them through matrix
for( i in 1:ncol(tmp_sub_means))
{
a <- pc %*% tmp_sub_means[ , i ]
pts_proj <- cbind(pts_proj, a)
}
# transpose matrix to long format
pts_proj <- t(pts_proj)
df <- as.data.frame( pts_proj)
circle fitting
sample <- pts_proj[ sample( nrow(pts_proj), size = 3, replace = F), ]
line_1_y_diff = sample[2, 2] - sample[1, 2]
line_1_x_diff = sample[2, 1] - sample[1, 1]
line_2_y_diff = sample[3, 2] - sample[2, 2]
line_2_x_diff = sample[3, 1] - sample[2, 1]
line_1_slope = line_1_y_diff/line_1_x_diff
line_2_slope = line_2_y_diff/line_2_x_diff
center_x = ( line_1_slope * line_2_slope *
( sample[1, 2] - sample[3, 2]) +
line_2_slope *
( sample[1, 1] + sample[2, 1]) -
line_1_slope *
( sample[2, 1] + sample[3, 1])) /
(2* (line_2_slope-line_1_slope) )
center_y = -1 * (center_x - ( sample[1, 1] + sample[2, 1])/2) /
line_1_slope +
(sample[1, 2]+sample[2, 2])/2
r = sqrt( (center_x - sample[1, 1])^2 + (center_y - sample[1, 2])^2)
center_x
center_y
r
ggplot( df, aes( x = pc2, y = pc3))+
geom_point()+
ggtitle("tree points projected")+
geom_circle( aes( x0 = center_y, y0 = center_x, r = r, color = "red"), inherit.aes = F )+
coord_fixed()+
theme( legend.position = "none")
sample <- pts_proj[ sample( nrow(pts_proj), size = 3, replace = F), ]
t0 <- Sys.time()
fit <- rcpp_circle_fit(sample)
Sys.time() - t0
fit
ggplot( df, aes( x = pc2, y = pc3))+
geom_point()+
ggtitle("tree points projected")+
geom_circle( aes( x0 = fit[1], y0 = fit[2], r = fit[3], color = "red"), inherit.aes = F )+
coord_fixed()+
theme( legend.position = "none")
tommywhitney/mobile_lidar documentation built on May 12, 2022, 7:38 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.
library(tidyverse) library(scatterplot3d) library(rgl) library(rospca) library(geometry) library(spanner) library(dbscan) library(ggforce) library(cowplot) library(Rcpp) library(mobileLidar) library(microbenchmark) # las <- readTLSLAS(file.choose()) # las <- readTLSLAS("C:/Users/matt/Desktop/DensePatchA.laz") # Rudd tank burn only north # las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_BO_N_zebcam_ALSregistered.laz", filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431804.383088774 3899195.49902686 30 (x y radius)") # rudd tank burn only south east # las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_BO_SE_zebcam_ALSregistered.laz", filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431838.287584286 3899147.61210088 30 (x y radius)") # # rudd tank burn only south west las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_BO_SW_zebcam_ALSregistered.laz", select = "xyz", filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431790.667864794 3899145.91081599 30 (x y radius)") # control north # las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_CO_N_zebcam_ALSregistered.laz", # select = "xyz", # filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431554.241897831 3899192.99624304 30 (x y radius)") # # control south east # las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_CO_SE_zebcam_ALSregistered.laz", # select = "xyz", # filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431595.891142449 3899156.16591304 30 (x y radius)") # control south west # las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_CO_SW_zebcam_ALSregistered.laz", # select = "xyz", # filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431545.984681478 3899149.36050936 30 (x y radius)") # treat and burn north east # las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_TB_NE_zebcam_ALSregistered.laz", # select = "xyz", # filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431940.528400044 3899570.10647836 30 (x y radius)") # treat and burn north west # las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_TB_NW_zebcam_ALSregistered.laz", # select = "xyz", # filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431883.255561657 3899552.23327632 30 (x y radius)") # # treat and burn south # las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_TB_S_zebcam_ALSregistered.laz", # select = "xyz", # filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431922.524931294 3899521.26763958 30 (x y radius)") # # treat north # las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_TO_N_zebcam_ALSregistered.laz", # select = "xyz", # filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431549.858895073 3899580.94101544 30 (x y radius)") # # treat south east # las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_TO_SE_zebcam_ALSregistered.laz", # select = "xyz", # filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431572.236974267 3899519.5420316 30 (x y radius)") # # treat south west # las <- readTLSLAS("C:/Users/matt/Desktop/Rudd_MLS_GeoSlam_ALSRegistered/Rudd_TO_SW_zebcam_ALSregistered.laz", # select = "xyz", # filter = "-keep_first -thin_with_voxel 0.05 -keep_circle 431506.243392164 3899537.98026232 30 (x y radius)") las <- classify_ground(las, csf(sloop_smooth = FALSE, class_threshold = 0.5, cloth_resolution = 0.5, rigidness = 1L, iterations = 500L, time_step = 0.65)) las <- normalize_height(las, tin()) las <- classify_noise(las, ivf(6,1)) las <- filter_poi(las, Classification != LASNOISE) plot(las, axis = T) las_slice <- filter_poi(las, Z>=0.62, Z<=2.12) eigens <- spanner::eigen_metrics(las_slice, radius = 0.33, ncpu = 8) pt_den <- spanner:::C_count_in_sphere(las_slice, radius = 0.33, ncpu = 8) las_slice@data<-cbind(las_slice@data, eigens) las_slice@data<-cbind(las_slice@data, pt_den) # adjust this Z slice if the trees width is close to the height of the slice las_slice <- filter_poi(las_slice, Z>=0.62 ,Z<= 2.12) clust <- dbscan::dbscan(las_slice@data[,c("X","Y","Z", "eSum","Verticality")], eps = 0.25, minPts = 100) las_slice@data$treeID<-clust$cluster plot(las_slice, color="treeID", axis = T)
fit_df <- las_slice_circle_fitting( las_slice, 80, 0.05, 0.6)
fit_df <- fit_df[fit_df$Radius <= 2, ] offsets <- plot(las_slice, color="treeID", axis = T) spheres3d( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = 1.37, r = fit_df[,3], alpha = .7) rgl.texts( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = fit_df[,3]+1.37 +.1, text = paste( "x: ",fit_df[,1])) rgl.texts( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = fit_df[,3]+1.37 +.2, text = paste( "y: ", fit_df[,2] )) rgl.texts( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = fit_df[,3]+1.37 +.3, text = paste( "diameter: ", fit_df[,3]*2)) rgl.texts( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = fit_df[,3]+1.37 +.4, text = paste( "mean squared error: ", fit_df[,4] )) rgl.texts( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = fit_df[,3]+1.37 +.5, text = paste( "inclusion: ", fit_df[,5] )) rgl.texts( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = fit_df[,3]+1.37 +.6, text = paste( "tree id: ", fit_df[,6] )) rgl.texts( x = fit_df[,1]-offsets[1], y = fit_df[,2]-offsets[2], z = fit_df[,3]+1.37 +.7, text = paste( "lean: ", fit_df[,7] ))
test_tree = 1 tree_slice <- filter_poi(las_slice, treeID == test_tree) # call las_slice_circle_fitting on tree slice to get a new fit with given parameters new_fit <- las_slice_circle_fitting( tree_slice, 80, 0.01, 0.6 ) fit_df <- rbind(fit_df, new_fit) plot(tree_slice, color="treeID", axis = T) # assign all points to a matrix M <- cbind( tree_slice$X, tree_slice$Y, tree_slice$Z, tree_slice$treeID) # select a tree from the matrix by id # tree_pts <- M[ M[ ,4] == i, ] # assign tmp matrix the points for current tree id tmp <- M[ , 1:3] # preform robust pca resR <- robpca(tmp, k=3, ndir = 5000) # filter tmp to around dbh # tmp <- tmp[ tmp[ ,3] >= 1.27, ] # tmp <- tmp[ tmp[ ,3] <= 1.47, ] # use robust pca to reduce dimensions and fit points to a plane orthogonal to the tree lean # get mean and principal components of points # note here pc1 will be the greatest value of the principal componets # with pc1 being the trees vertical componet, and pc2 and pc3 being the horizontal components # if the tree has a greater width than the height of the slice taken, what we assume to be pc1 as # the height will be one of the horozontal componets center <- resR$center pc1 <- resR$loadings[ , 1] pc2 <- resR$loadings[ , 2] pc3 <- resR$loadings[ , 3] # bind the lower two principal components into a matrix pc <- t(cbind( pc3, pc2)) # subtract means from each point tmp_sub_means <- t( sweep( tmp, 2, center )) # init empty matrix pts_proj <- matrix( 0,nrow = 2) # loop through all points passing them through matrix for( i in 1:ncol(tmp_sub_means)) { a <- pc %*% tmp_sub_means[ , i ] pts_proj <- cbind(pts_proj, a) } # for some reason when removing these na values the loop does not end # remove first row of matrix, contains na values # pts_proj <-pts_proj[ -1,] # transpose matrix to long format pts_proj <- t(pts_proj) pts_proj <- pts_proj[ -c(1),] df <- as.data.frame( pts_proj) best_fit <- ransac_circle_fit( pts_proj, 10000, 0.1, 0.8) format(best_fit, scientific = F) center xyz_center <- center + best_fit[1]%*%pc3 + best_fit[2]%*%pc2 plot(tree_slice, color="treeID") # spheres3d(x = xyz_center[1], y = xyz_center[2], z = 1.37, radius = best_fit[3]) ggplot( df, aes( x = pc3, y = pc2))+ geom_point()+ ggtitle("tree points projected")+ geom_circle( aes( x0 = best_fit[1], y0 = best_fit[2], r = best_fit[3], color = "red"), inherit.aes = F )+ coord_fixed()+ theme( legend.position = "none")+ geom_point(aes(x=best_fit[1], y=best_fit[2]), colour="blue")+ coord_fixed()
offsets <- plot(tree_slice, color="treeID") spheres3d(x = xyz_center[,1]-offsets[1], y = xyz_center[,2]-offsets[2], z = 1.37, r = fit_df[,3], alpha = .7) # rgl.texts( x = xyz_center[,1]-offsets[1], y = xyz_center[,2]-offsets[2], z = 2.3, text = paste( "diameter: ", fit_df[,3]*2)) # rgl.texts( x = xyz_center[,1]-offsets[1], y = xyz_center[,2]-offsets[2], z = 2.15, text = paste( "x: ",fit_df[,1])) # rgl.texts( x = xyz_center[,1]-offsets[1], y = xyz_center[,2]-offsets[2], z = 2, text = paste( "y: ", fit_df[,2] ))
open3d() theta <- seq(0, 2*pi, len = 25) knot <- cylinder3d( center = cbind( sin(theta) + 2*sin(2*theta), 2*sin(3*theta), cos(theta) - 2*cos(2*theta)), e1 = cbind( cos(theta) + 4*cos(2*theta), 6*cos(3*theta), sin(theta) + 4*sin(2*theta)), radius = 0.8, closed = TRUE)
# assign all points to a matrix M <- cbind( tree_slice$X, tree_slice$Y, tree_slice$Z, tree_slice$treeID) # assign tmp matrix the points for current tree id tmp <- M[ , 1:3] # preform robust pca resR <- robpca(tmp, k=3, ndir = 5000) # filter tmp to around dbh # tmp <- tmp[ tmp[ ,3] >= 1.27, ] # tmp <- tmp[ tmp[ ,3] <= 1.47, ] # use robust pca to reduce dimensions and fit points to a plane orthogonal to the tree lean # get mean and principal components of points # note here pc1 will be the greatest value of the principal componets # with pc1 being the trees vertical componet, and pc2 and pc3 being the horizontal components # if the tree has a greater width than the height of the slice taken, what we assume to be pc1 as # the height will be one of the horozontal componets center <- resR$center pc1 <- resR$loadings[ , 1] pc2 <- resR$loadings[ , 2] pc3 <- resR$loadings[ , 3] # bind the lower two principal components into a matrix pc <- t(cbind( pc3, pc2)) # subtract means from each point tmp_sub_means <- t( sweep( tmp, 2, center )) # init empty matrix pts_proj <- matrix( 0,nrow = 2) # loop through all points passing them through matrix for( i in 1:ncol(tmp_sub_means)) { a <- pc %*% tmp_sub_means[ , i ] pts_proj <- cbind(pts_proj, a) } # for some reason when removing these na values the loop does not end # remove first row of matrix, contains na values # pts_proj <-pts_proj[ -1,] # transpose matrix to long format pts_proj <- t(pts_proj) pts_proj <- pts_proj[ -c(1),] df <- as.data.frame( pts_proj) ggplot( df, aes( x = pc3, y = pc2))+ geom_point()+ ggtitle("tree points projected")+ geom_circle( aes( x0 = fit_df[test_tree,1], y0 = fit_df[test_tree,2], r = fit_df[test_tree,3], color = "red"), inherit.aes = F )+ coord_fixed()+ theme( legend.position = "none")+ coord_fixed() fit_df[test_tree,]
# assign all points to a matrix M <- cbind( las_slice$X, las_slice$Y, las_slice$Z, las_slice$treeID) # select a tree from the matrix by id tree_pts <- M[ M[ ,4] == i, ] # assign tmp matrix the points for current tree id tmp <- tree_pts[ , 1:3] # preform robust pca resR <- robpca(tmp, k=3, ndir = 5000) # filter tmp to around dbh # tmp <- tmp[ tmp[ ,3] >= 1.27, ] # tmp <- tmp[ tmp[ ,3] <= 1.47, ] # use robust pca to reduce dimensions and fit points to a plane orthogonal to the tree lean # get mean and principal components of points # note here pc1 will be the greatest value of the principal componets # with pc1 being the trees vertical componet, and pc2 and pc3 being the horizontal components # if the tree has a greater width than the height of the slice taken, what we assume to be pc1 as # the height will be one of the horozontal componets center <- resR$center pc1 <- resR$loadings[ , 1] pc2 <- resR$loadings[ , 2] pc3 <- resR$loadings[ , 3] # bind the lower two principal components into a matrix pc <- t(cbind( pc3, pc2)) # subtract means from each point tmp_sub_means <- t( sweep( tmp, 2, center )) # init empty matrix pts_proj <- matrix( 0,nrow = 2) # loop through all points passing them through matrix for( i in 1:ncol(tmp_sub_means)) { a <- pc %*% tmp_sub_means[ , i ] pts_proj <- cbind(pts_proj, a) } # for some reason when removing these na values the loop does not end # remove first row of matrix, contains na values # pts_proj <-pts_proj[ -1,] # transpose matrix to long format pts_proj <- t(pts_proj) pts_proj <- pts_proj[ -c(1),] df <- as.data.frame( pts_proj) best_fit <- ransac_circle_fit( pts_proj, 1000, 0.01, 0.8) format(best_fit, scientific = F) best_fit[3]*2 geom_point()+ ggtitle("tree points projected")+ geom_circle( aes( x0 = best_fit[1], y0 = best_fit[2], r = best_fit[3], color = "red"), inherit.aes = F )+ coord_fixed()+ theme( legend.position = "none")+ geom_point(aes(x=best_fit[1], y=best_fit[2]), colour="blue")+ coord_fixed()
rgl.open() rgl.points(tree_pts[ , 1 ], tree_pts[ , 2], tree_pts[ , 3], color ="white", axis = T)
# select a tree from the matrix by id tree_pts <- M[ M[ ,4] == 1, ] # assign tmp matrix the points for current tree id tmp <- tree_pts[ , 1:3] # preform robust pca resR <- robpca(tmp, k=3, ndir = 5000) # use robust pca to reduce dimensions and fit points to a plane orthogonal to the tree lean # get mean and principal components of points # note here pc1 will be the greatest value of the principal componets # with pc1 being the trees vertical componet, and pc2 and pc3 being the horizontal components # if the tree has a greater width than the height of the slice taken, what we assume to be pc1 as # the height will be one of the horozontal componets center <- resR$center pc1 <- resR$loadings[ , 1] pc2 <- resR$loadings[ , 2] pc3 <- resR$loadings[ , 3] # bind the lower two principal components into a matrix pc <- t(cbind( pc3, pc2))(pts_proj) # subtract means from each point tmp_sub_means <- t( sweep( tmp, 2, center )) # init empty matrix pts_proj <- matrix( nrow = 2) # loop through all points passing them through matrix for( i in 1:ncol(tmp_sub_means)) { a <- pc %*% tmp_sub_means[ , i ] pts_proj <- cbind(pts_proj, a) } # transpose matrix to long format pts_proj <- t df <- as.data.frame( pts_proj) # }
best_fit <- ransac_circle_fit( pts_proj, 10000, 0.001, 0.9) best_fit
circle fitting
ggplot( df, aes( x = pc3, y = pc2))+ geom_point()+ ggtitle("tree points projected")+ geom_circle( aes( x0 = best_fit[1], y0 = best_fit[2], r = best_fit[3], color = "red"), inherit.aes = F )+ coord_fixed()+ theme( legend.position = "none") best_fit
best_fit[3]*2
fit at 20cm around dbh
# assign all points to a matrix M <- cbind( las_slice$X, las_slice$Y, las_slice$Z, las_slice$treeID) # select a tree from the matrix by id tree_pts <- M[M[ ,4] == 1,] # assign tmp matrix the points for current tree id tmp <- tree_pts[ , 1:3] # preform robust pca resR <- robpca(tmp, k=3, ndir = 5000) # select 10 cm of points around dbh tmp <- tmp[ tmp[ , 3] >= 1.27 & tmp[ , 3] <= 1.47, ] # use robust pca to reduce dimensions and fit points to a plane orthogonal to the tree lean # get mean and principal components of points # note here pc1 will be the greatest value of the principal componets # with pc1 being the trees vertical componet, and pc2 and pc3 being the horizontal componets # if the tree has a greater width than the height of the slice taken what we assume to be pc1 as # the height will be one of the horozontal componets center <- resR$center pc1 <- resR$loadings[ , 1] pc2 <- resR$loadings[ , 2] pc3 <- resR$loadings[ , 3] # bind the lower two principal components into a matrix pc <- t(cbind( pc3, pc2)) # subtract means from each point tmp_sub_means <- t( sweep( tmp, 2, center )) # init empty matrix pts_proj <- matrix( nrow = 2) # loop through all points passing them through matrix for( i in 1:ncol(tmp_sub_means)) { a <- pc %*% tmp_sub_means[ , i ] pts_proj <- cbind(pts_proj, a) } # transpose matrix to long format pts_proj <- t(pts_proj) df <- as.data.frame( pts_proj)
best_fit <- ransac_circle_fit( pts_proj, 10000, 0.001, 0.9) best_fit
ggplot( df, aes( x = pc2, y = pc3))+ geom_point()+ ggtitle("tree points projected")+ geom_circle( aes( x0 = best_fit[2], y0 = best_fit[1], r = best_fit[3], color = "red"), inherit.aes = F )+ coord_fixed()+ theme( legend.position = "none")
circle fitting
sample <- pts_proj[ sample( nrow(pts_proj), size = 3, replace = F), ]
line_1_y_diff = sample[2, 2] - sample[1, 2] line_1_x_diff = sample[2, 1] - sample[1, 1] line_2_y_diff = sample[3, 2] - sample[2, 2] line_2_x_diff = sample[3, 1] - sample[2, 1] line_1_slope = line_1_y_diff/line_1_x_diff line_2_slope = line_2_y_diff/line_2_x_diff center_x = ( line_1_slope * line_2_slope * ( sample[1, 2] - sample[3, 2]) + line_2_slope * ( sample[1, 1] + sample[2, 1]) - line_1_slope * ( sample[2, 1] + sample[3, 1])) / (2* (line_2_slope-line_1_slope) ) center_y = -1 * (center_x - ( sample[1, 1] + sample[2, 1])/2) / line_1_slope + (sample[1, 2]+sample[2, 2])/2 r = sqrt( (center_x - sample[1, 1])^2 + (center_y - sample[1, 2])^2) center_x center_y r
view single iteration circle fit
ggplot( df, aes( x = pc2, y = pc3))+ geom_point()+ ggtitle("tree points projected")+ geom_circle( aes( x0 = center_y, y0 = center_x, r = r, color = "red"), inherit.aes = F )+ coord_fixed()+ theme( legend.position = "none")
# assign all points to a matrix M <- cbind( las_slice$X, las_slice$Y, las_slice$Z, las_slice$treeID) # select a tree from the matrix by id tree_pts <- M[M[ ,4] == 2,] # assign tmp matrix the points for current tree id tmp <- tree_pts[ , 1:3] # preform robust pca resR <- robpca(tmp, k=3, ndir = 5000) # select 10 cm of points around dbh tmp <- tmp[ tmp[ , 3] >= 1.27 & tmp[ , 3] <= 1.47, ] # use robust pca to reduce dimensions and fit points to a plane orthogonal to the tree lean # get mean and principal components of points # note here pc1 will be the greatest value of the principal componets # with pc1 being the trees vertical componet, and pc2 and pc3 being the horizontal componets # if the tree has a greater width than the height of the slice taken what we assume to be pc1 as # the height will be one of the horozontal componets center <- resR$center pc1 <- resR$loadings[ , 1] pc2 <- resR$loadings[ , 2] pc3 <- resR$loadings[ , 3] # bind the lower two principal components into a matrix pc <- t(cbind( pc3, pc2)) # subtract means from each point tmp_sub_means <- t( sweep( tmp, 2, center )) # init empty matrix pts_proj <- matrix( nrow = 2) # loop through all points passing them through matrix for( i in 1:ncol(tmp_sub_means)) { a <- pc %*% tmp_sub_means[ , i ] pts_proj <- cbind(pts_proj, a) } # transpose matrix to long format pts_proj <- t(pts_proj) df <- as.data.frame( pts_proj)
circle fitting
sample <- pts_proj[ sample( nrow(pts_proj), size = 3, replace = F), ]
line_1_y_diff = sample[2, 2] - sample[1, 2] line_1_x_diff = sample[2, 1] - sample[1, 1] line_2_y_diff = sample[3, 2] - sample[2, 2] line_2_x_diff = sample[3, 1] - sample[2, 1] line_1_slope = line_1_y_diff/line_1_x_diff line_2_slope = line_2_y_diff/line_2_x_diff center_x = ( line_1_slope * line_2_slope * ( sample[1, 2] - sample[3, 2]) + line_2_slope * ( sample[1, 1] + sample[2, 1]) - line_1_slope * ( sample[2, 1] + sample[3, 1])) / (2* (line_2_slope-line_1_slope) ) center_y = -1 * (center_x - ( sample[1, 1] + sample[2, 1])/2) / line_1_slope + (sample[1, 2]+sample[2, 2])/2 r = sqrt( (center_x - sample[1, 1])^2 + (center_y - sample[1, 2])^2) center_x center_y r
ggplot( df, aes( x = pc2, y = pc3))+ geom_point()+ ggtitle("tree points projected")+ geom_circle( aes( x0 = center_y, y0 = center_x, r = r, color = "red"), inherit.aes = F )+ coord_fixed()+ theme( legend.position = "none")
sample <- pts_proj[ sample( nrow(pts_proj), size = 3, replace = F), ] t0 <- Sys.time() fit <- rcpp_circle_fit(sample) Sys.time() - t0 fit ggplot( df, aes( x = pc2, y = pc3))+ geom_point()+ ggtitle("tree points projected")+ geom_circle( aes( x0 = fit[1], y0 = fit[2], r = fit[3], color = "red"), inherit.aes = F )+ coord_fixed()+ theme( legend.position = "none")
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.