diff --git a/CodeR/01_test_data_structure.r b/CodeR/01_Compute_angles_square.r similarity index 100% rename from CodeR/01_test_data_structure.r rename to CodeR/01_Compute_angles_square.r diff --git a/CodeR/02_Compute_angles_cone.r b/CodeR/02_Compute_angles_cone.r new file mode 100644 index 0000000000000000000000000000000000000000..754d4f0ab512d44a23a4b622b8317daa4c1fddee --- /dev/null +++ b/CodeR/02_Compute_angles_cone.r @@ -0,0 +1,122 @@ +rm(list=ls()) +gc() + +## system("gio mount -u smb://int.nibio.no/Databank") +## system(paste0("gio mount smb://int.nibio.no/Databank ", +## "< /home/nibio/.craden")) +## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +## Move the scan to have a local copy. +## Origin <- paste0("/run/user/1003/gvfs", +## "/smb-share:server=int.nibio.no,", +## "share=databank/Prosjekter", +## "/51160_IMPRINT/Forest_projections", +## "/DataR/to_map/merged_22_classified.las") +Dest <- "/home/nibio/DATA" +## cmd <- paste("cp", Origin, Dest) +## system(cmd) + +## REMOVE FILE FROM ORIGIN! + +## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +rm(list=ls()) + +library("rlas") +library("rgl") +library("stringr") +library("VoxR") +install.packages("VoxR") + +install.packages("Rfast") + +library("data.table") + +load("/home/pepito/Documents/holger_tower_pov/DataR/pan3d.RData") + +## this function scale points coordinate +## between 0 and (max-min) values +Scale <- function(Vector){ + ((Vector - min(Vector))/ + (max(Vector)-min(Vector)))* + (max(Vector)-min(Vector))} + +## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +## plot data - full path +DataPaths <- list.files(paste0("/home/pepito/Documents", + "/holger_tower_pov/DataR"), + pattern = "plot_", + ## pattern = "Plot_", + full.names = T) + +Data <- read.las(DataPaths) +Data <- Data[sample(1:nrow(Data), nrow(Data)*0.20), ] + +Data[, X := Scale(Data$X)] +Data[, Y := Scale(Data$Y)] +Data[, Z := Scale(Data$Z)] + +Data[, X := X-mean(X)] +Data[, Y := Y-mean(Y)] +Data <- Data[, 1:3] + + +plot3d(Data, aspect = "iso", size = 0.5) +pan3d(3) +## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +## Create an artificial tower. Just the corners +## CornA <- data.table(X = 3, Y = 3, Z = 40) +CornA <- Data[identify3d(Data, n = 1), ] +spheres3d(CornA, radius = 1, col = "yellow") + +## Left click on a ground point in the point cloud +## Right click to scape +CornB <- Data[identify3d(Data, n = 1), ] +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), + by = seq_len(nrow(Data))] +setnames(Data, "Dista", "CO") + +## plot3d(Data, aspect = "iso", size = 0.5, col = "grey21") +## plot3d(Data[CO < 2 , ], aspect = "iso", size = 4, col = "red", add = T) +## spheres3d(CornA, radius = 1, col = "yellow") +## spheres3d(CornB, radius = 1, col = "red") +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), + by = seq_len(nrow(Data))] + +## plot3d(Data2, aspect = "iso", size = 4, col = "red", add = T) +## spheres3d(CornA, radius = 1, col = "yellow") +## spheres3d(CornB, radius = 1, col = "red") +Data[, Angle := CO/H] +Data[, Angle := (asin(Angle))*180/pi] + +plot3d(Data[Angle < 7, ], aspect = "iso", size = 4, col = "red", add = T) +pan3d(3) diff --git a/CodeR/02_Compute_angles_cone.r~ b/CodeR/02_Compute_angles_cone.r~ new file mode 100644 index 0000000000000000000000000000000000000000..56a7d3fef265ce5a31415d68ab737368d92069bb --- /dev/null +++ b/CodeR/02_Compute_angles_cone.r~ @@ -0,0 +1,118 @@ +rm(list=ls()) +gc() + +## system("gio mount -u smb://int.nibio.no/Databank") +## system(paste0("gio mount smb://int.nibio.no/Databank ", +## "< /home/nibio/.craden")) +## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +## Move the scan to have a local copy. +## Origin <- paste0("/run/user/1003/gvfs", +## "/smb-share:server=int.nibio.no,", +## "share=databank/Prosjekter", +## "/51160_IMPRINT/Forest_projections", +## "/DataR/to_map/merged_22_classified.las") +Dest <- "/home/nibio/DATA" +## cmd <- paste("cp", Origin, Dest) +## system(cmd) + +## REMOVE FILE FROM ORIGIN! + +## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +rm(list=ls()) + +library("rlas") +library("rgl") +library("stringr") +library("VoxR") +library("data.table") + +load("/home/pepito/Documents/holger_tower_pov/DataR/pan3d.RData") + +## this function scale points coordinate +## between 0 and (max-min) values +Scale <- function(Vector){ + ((Vector - min(Vector))/ + (max(Vector)-min(Vector)))* + (max(Vector)-min(Vector))} + +## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +## plot data - full path +DataPaths <- list.files(paste0("/home/pepito/Documents", + "/holger_tower_pov/DataR"), + pattern = "plot_", + ## pattern = "Plot_", + full.names = T) + +Data <- read.las(DataPaths) +Data <- Data[sample(1:nrow(Data), nrow(Data)*0.20), ] + +Data[, X := Scale(Data$X)] +Data[, Y := Scale(Data$Y)] +Data[, Z := Scale(Data$Z)] + +Data[, X := X-mean(X)] +Data[, Y := Y-mean(Y)] +Data <- Data[, 1:3] + + +plot3d(Data, aspect = "iso", size = 0.5) +pan3d(3) +## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +## Create an artificial tower. Just the corners +## CornA <- data.table(X = 3, Y = 3, Z = 40) +CornA <- Data[identify3d(Data, n = 1), ] +spheres3d(CornA, radius = 1, col = "yellow") + +## Left click on a ground point in the point cloud +## Right click to scape +CornB <- Data[identify3d(Data, n = 1), ] +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), + by = seq_len(nrow(Data))] +setnames(Data, "Dista", "CO") + +## plot3d(Data, aspect = "iso", size = 0.5, col = "grey21") +## plot3d(Data[CO < 2 , ], aspect = "iso", size = 4, col = "red", add = T) +## spheres3d(CornA, radius = 1, col = "yellow") +## spheres3d(CornB, radius = 1, col = "red") +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), + by = seq_len(nrow(Data))] + +## plot3d(Data2, aspect = "iso", size = 4, col = "red", add = T) +## spheres3d(CornA, radius = 1, col = "yellow") +## spheres3d(CornB, radius = 1, col = "red") +Data[, Angle := CO/H] +Data[, Angle := (asin(Angle))*180/pi] + +plot3d(Data[Angle < 7, ], aspect = "iso", size = 4, col = "red", add = T) +pan3d(3)