diff --git a/CodeR/02_Compute_angles_cone.r b/CodeR/02_Compute_angles_cone.r index 555097c0a3dad02ec40474552119cfc31a65aee1..19ea241eff0cb6e4e132fa03f0f5940da0fdd706 100644 --- a/CodeR/02_Compute_angles_cone.r +++ b/CodeR/02_Compute_angles_cone.r @@ -16,6 +16,7 @@ gc() ## system(cmd) ## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +## packages and functions rm(list=ls()) library("rlas") @@ -34,24 +35,54 @@ Scale <- function(Vector){ (max(Vector)-min(Vector)))* (max(Vector)-min(Vector))} +## Distances +cross3d_prod <- function(v1,v2){ + v3 <- vector() + v3[1] <- v1[2]*v2[3]-v1[3]*v2[2] + v3[2] <- v1[3]*v2[1]-v1[1]*v2[3] + v3[3] <- v1[1]*v2[2]-v1[2]*v2[1] + return(v3) +} + +dist3d <- function(xa, ya, za, xb, yb, zb, xc, yc, zc) { + a <- c(xa, ya, za) + b <- c(xb, yb, zb, xc) + c <- c(xc, yc, zc) + + v1 <- b - c + v2 <- a - b + v3 <- cross3d_prod(v1,v2) + area <- sqrt(sum(v3*v3))/2 + d <- 2*area/sqrt(sum(v1*v1)) + return(d) +} + +dist3d2 <- function(xa, ya, za, xb, yb, zb) { + d <- sqrt(((xa - xb)^2) + ((ya - yb)^2) + ((za - zb)^2)) + return(d)} + ## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ## DataR DataR <- paste0("/home/nibio/REPOS/holger_tower_pov/DataR") ## DataR <- paste0("/home/pepito/Documents/holger_tower_pov/DataR") + + + ## plot data - full path DataPaths <- list.files(DataR, ## pattern = "plot_", ## pattern = "Plot_", ## pattern = "merg", - pattern = ".rds", + ## pattern = ".rds", + pattern = "Tower", full.names = T) -## Data <- read.las(DataPaths) -## Data[, 1:3] -Data <- fst::read_fst(DataPaths, as.data.table = T) +Data <- read.las(DataPaths) +Data[, 1:3] +## Data <- fst::read_fst(DataPaths, as.data.table = T) -Data <- Data[sample(1:nrow(Data), nrow(Data)*0.02), ] +Data <- Data[sample(1:nrow(Data), nrow(Data)*0.10), ] Data[, X := Scale(Data$X)] Data[, Y := Scale(Data$Y)] Data[, Z := Scale(Data$Z)] @@ -59,9 +90,19 @@ Data[, Z := Scale(Data$Z)] Data[, X := X-mean(X)] Data[, Y := Y-mean(Y)] + +HGSTPoint <- Data[max(Z), ] + +Data[, Dist := dist3d2(X, Y, 0, + HGSTPoint$X, HGSTPoint$Y, 0), + by = seq_len(nrow(Data))] + + ## open3d() -Data <- Data[Classification != 2,] -plot3d(Data, aspect = "iso", size = 0.5) +Data <- Data[Dist < 10,] +plot3d(Data, aspect = "iso", size = 3) +spheres3d(HGSTPoint, radius = 1, col = "yellow") + pan3d(3) ## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ## Create an artificial tower. Just the corners @@ -75,27 +116,6 @@ spheres3d(CornA, radius = 1, col = "yellow") spheres3d(CornB, radius = 1, col = "red") ## rgl.pop() -cross3d_prod <- function(v1,v2){ - v3 <- vector() - v3[1] <- v1[2]*v2[3]-v1[3]*v2[2] - v3[2] <- v1[3]*v2[1]-v1[1]*v2[3] - v3[3] <- v1[1]*v2[2]-v1[2]*v2[1] - return(v3) -} - -dist3d <- function(xa, ya, za, xb, yb, zb, xc, yc, zc) { - a <- c(xa, ya, za) - b <- c(xb, yb, zb, xc) - c <- c(xc, yc, zc) - - v1 <- b - c - v2 <- a - b - v3 <- cross3d_prod(v1,v2) - area <- sqrt(sum(v3*v3))/2 - d <- 2*area/sqrt(sum(v1*v1)) - return(d) -} - Data[, Dista := dist3d(X, Y, Z, CornA$X, CornA$Y, CornA$Z, CornB$X, CornB$Y, CornB$Z), @@ -107,9 +127,6 @@ setnames(Data, "Dista", "CO") ## spheres3d(CornA, radius = 1, col = "yellow") ## spheres3d(CornB, radius = 1, col = "red") ## pan3d(3) -dist3d2 <- function(xa, ya, za, xb, yb, zb) { - d <- sqrt(((xa - xb)^2) + ((ya - yb)^2) + ((za - zb)^2)) - return(d)} Data[, H := dist3d2(X, Y, Z, CornA$X, CornA$Y, CornA$Z),