diff --git a/CodeR/03_Compute_angles_cone.r b/CodeR/03_Compute_angles_cone.r
index eb0af4d713cfba0199cef2c1a29574101f4b0ada..c312c9b821d7bcd268da3617db3254c4b77d447e 100644
--- a/CodeR/03_Compute_angles_cone.r
+++ b/CodeR/03_Compute_angles_cone.r
@@ -18,13 +18,18 @@ library("data.table")
## DataR <- paste0("/home/nica/Documents/Holger_tower_pov/DataR")
## CodeR <- paste0("/home/nica/Documents/Holger_tower_pov/CodeR")
-DataR <- paste0("/home/nibio/REPOS/holger_tower_pov/DataR")
-CodeR <- paste0("/home/nibio/REPOS/holger_tower_pov/CodeR")
+## DataR <- paste0("/home/nibio/REPOS/holger_tower_pov/DataR")
+## CodeR <- paste0("/home/nibio/REPOS/holger_tower_pov/CodeR")
+
+CodeR <- "/home/pepito/Documents/holger_tower_pov/CodeR/"
+DataR <- "/home/pepito/Desktop/Temp_Holger_Data/"
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## Functions
source(paste0(CodeR, "/Functions/Functions.r"))
+objects()
+dist3d
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## plot data - full path
DataPaths <- list.files(DataR,
@@ -35,159 +40,145 @@ DataPaths <- list.files(DataR,
## pattern = "Tower",
full.names = T)
-Data <- read.las(DataPaths[1])
+## Data <- read.las(DataPaths[1])
+Data <- fst::read_fst(DataPaths[1], as.data.table = T)
Data <- Data[, 1:3]
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## Hi Nico,
+## Assuming the coordinates in the PointCloud are UTM,
+## the preferred location of the SIF sensor would be at
-# Create a simple 3D arrow function
-drawArrow <- function(start, end, col = "red", lwd = 5, arrowSize = 0.1) {
- segments3d(rbind(start, end), col = col, lwd = lwd)
-
- # Compute the direction vector
- dirVector <- end - start
-
- # Normalize the direction vector
- dirVector <- dirVector / sqrt(sum(dirVector^2))
-
- # Create the arrow head (as three segments emanating from the arrow tip)
- segments3d(rbind(end, end - arrowSize * dirVector + c( arrowSize, 0, 0)), col = col, lwd = lwd)
- segments3d(rbind(end, end - arrowSize * dirVector + c(-arrowSize, 0, 0)), col = col, lwd = lwd)
- segments3d(rbind(end, end - arrowSize * dirVector + c(0, -arrowSize, 0)), col = col, lwd = lwd)
-}
-
-# Example data
-set.seed(123)
-n <- 1000
-x <- runif(n, 500000, 505000) # UTM x-coordinates
-y <- runif(n, 5000000, 5005000) # UTM y-coordinates
-z <- runif(n, 0, 100) # just some random z values
-
-# Plot the point cloud
-plot3d(x, y, z, type = "p", col = "blue")
-
-# Define the "Devise" point
-devise <- c(502500, 5002500, 50) # For instance
-spheres3d(devise, radius = 100)
-
-
-# Compute the arrow length based on the data range
-arrow_length <- max(c(diff(range(x)), diff(range(y)), diff(range(z)))) / 10
-
-# South (S)
-drawArrow(devise, devise + c(0, -arrow_length, 0), col="red")
-
-# West (W)
-drawArrow(devise, devise + c(-arrow_length, 0, 0), col="green")
-
-# North (N)
-drawArrow(devise, devise + c(0, arrow_length, 0), col="blue")
-
-# East (E)
-drawArrow(devise, devise + c(arrow_length, 0, 0), col="yellow")
-
-# Southwest (SW)
-drawArrow(devise, devise + c(-arrow_length, -arrow_length, 0), col="purple")
+## N 6694611.571 E 614675.0338 and height 315.8627 masl
-# Southeast (SE)
-drawArrow(devise, devise + c(arrow_length, -arrow_length, 0), col="cyan")
+## This is a point in the middle of the SW-SE edge of the tower
+## top, and the preferred viewing direction is more or less south
+## since there is a dense group of trees there; the crucial
+## question is the viewing angle.
-# Northwest (NW)
-drawArrow(devise, devise + c(-arrow_length, arrow_length, 0), col="magenta")
+## Cheers, Holger
+## Devise <- data.table(X = 614675.0338,
+## Y = 6694611.571,
+## Z = 315.8627)
-# Northeast (NE)
-drawArrow(devise, devise + c(arrow_length, arrow_length, 0), col="orange")
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+Data <- Data[sample(1:nrow(Data), nrow(Data)*0.10), ]
+gc()
+## Data <- rbind(Data, Devise)
+Data[, X := Scale(Data$X)]
+Data[, Y := Scale(Data$Y)]
+## Devise <- tail(Data, 1)
-getArrowTipCoords <- function(base, length, direction) {
- if (direction == "S") {
- return(base + c(0, -length, 0))
- } else if (direction == "W") {
- return(base + c(-length, 0, 0))
- } else if (direction == "N") {
- return(base + c(0, length, 0))
- } else if (direction == "E") {
- return(base + c(length, 0, 0))
- } else if (direction == "SW") {
- return(base + c(-length, -length, 0))
- } else if (direction == "SE") {
- return(base + c(length, -length, 0))
- } else if (direction == "NW") {
- return(base + c(-length, length, 0))
- } else if (direction == "NE") {
- return(base + c(length, length, 0))
- } else {
- stop("Invalid direction provided!")
- }
-}
+plot3d(Data, aspect = "iso", size = 0.2)
+pan3d(3)
+## spheres3d(Devise, radius = 0.5, col = "red")
-# Define the "Devise" point and arrow length as before
-devise <- c(502500, 5002500, 50) # For instance
-arrow_length <- max(c(diff(range(x)), diff(range(y)), diff(range(z)))) / 10
+## Devise location is not in the area mentioned in the email.
-# Get the coordinates for each direction
-s_tip <- getArrowTipCoords(devise, arrow_length, "S")
-w_tip <- getArrowTipCoords(devise, arrow_length, "W")
-n_tip <- getArrowTipCoords(devise, arrow_length, "N")
-e_tip <- getArrowTipCoords(devise, arrow_length, "E")
-sw_tip <- getArrowTipCoords(devise, arrow_length, "SW")
-se_tip <- getArrowTipCoords(devise, arrow_length, "SE")
-nw_tip <- getArrowTipCoords(devise, arrow_length, "NW")
-ne_tip <- getArrowTipCoords(devise, arrow_length, "NE")
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## Select the top of the tower
+## ids <- plot3d(Data, aspect = "iso", size = 4)
+## pan3d(3)
+## TopTowerIndex <- selectpoints3d(ids["data"], value = FALSE)
-# Print the coordinates
-print(list(S = s_tip))
+## TopTowerIndex[, 2]
+## plot3d(Data[-TopTowerIndex[, 2], ], aspect = "iso", size = 3)
+## pan3d(3)
+## plot3d(Data[TopTowerIndex[, 2], ], aspect = "iso",
+## size = 3, col = "red", add = T)
+## plot3d(Data[TopTowerIndex[, 2], ], aspect = "iso",
+## size = 5, col = "red")
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## TopTowerP <- Data[TopTowerIndex[, 2], ]
+## saveRDS(TopTowerP, paste0(DataR, "TopTowerP.rds"))
+## TopTowerP <- readRDS(paste0(DataR, "TopTowerP.rds"))
+## plot3d(TopTowerP, aspect = "iso",
+## size = 10, col = "red")
+## ToWerCenter <- TopTowerP[10, ]
+## spheres3d(ToWerCenter, radius = 0.1, col = "black")
+## TopTowerSW <- TopTowerP[306, ]
+## spheres3d(TopTowerSW, radius = 0.1, col = "green")
+## TopTowerSE <- TopTowerP[168, ]
+## spheres3d(TopTowerSE, radius = 0.1, col = "blue")
-## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
-## Hi Nico,
-## Assuming the coordinates in the PointCloud are UTM,
-## the preferred location of the SIF sensor would be at
-## N 6694611.571 E 614675.0338 and height 315.8627 masl
-## This is a point in the middle of the SW-SE edge of the tower
-## top, and the preferred viewing direction is more or less south
-## since there is a dense group of trees there; the crucial
-## question is the viewing angle.
+cloud <- data.frame(x = runif(10000, -20, 20), y = runif(10000, -20, 20), z = runif(10000, 0, 40))
-## Cheers, Holger
+# Define the TopTowerSW point
+TopTowerSW <- data.frame(x = 0.5, y = 0.5, x= 40)
-Devise <- data.table(X = 614675.0338,
- Y = 6694611.571,
- Z = 315.8627)
+# Function to compute the distance between points and the line
+distance_to_line <- function(cloud, TopTowerSW, theta_deg) {
+ # Convert angle to radians
+ theta_rad <- theta_deg * pi / 180
+
+ # Direction vector for the line
+ d <- c(0, -1, tan(theta_rad))
+
+ # Filter the cloud to retain only points to the south of TopTowerSW
+ cloud_south <- cloud[cloud$y < TopTowerSW$y, ]
+
+ # Convert both data.frames to matrices for easier computation
+ cloud_matrix <- as.matrix(cloud_south)
+ TopTowerSW_matrix <- matrix(rep(as.matrix(TopTowerSW), nrow(cloud_south)), ncol=3, byrow=TRUE)
+
+ # Compute the position vectors relative to TopTowerSW for all points
+ r <- cloud_matrix - TopTowerSW_matrix
+
+ # Compute the cross product between r and d for all points
+ cross_product <- cbind(r[,2]*d[3] - r[,3]*d[2], r[,3]*d[1] - r[,1]*d[3], r[,1]*d[2] - r[,2]*d[1])
+
+ # Compute the distance using the formula mentioned
+ distances <- sqrt(rowSums(cross_product^2)) / sqrt(sum(d^2))
+
+ return(distances)
+}
-## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
-Data <- Data[sample(1:nrow(Data), nrow(Data)*0.10), ]
+# Example usage: Compute the distance with a vertical angle of 45 degrees
+distances <- distance_to_line(cloud, TopTowerSW, -60)
-Data <- rbind(Data, Devise)
-Data[, X := Scale(Data$X)]
-Data[, Y := Scale(Data$Y)]
+plot3d(cloud, size = 5)
+plot3d(TopTowerSW, size = 10, add = T, col = "red")
-Data[, X := X-mean(X)]
-Data[, Y := Y-mean(Y)]
+Puntitos <- as.numeric(names(distances[distances < 4]))
+plot3d(cloud[Puntitos, ], size = 10, col = "blue", add = T)
-Devise <- tail(Data, 1)
+compute_endpoint <- function(TopTowerSW, cloud, theta_deg) {
+ # Direction vector for the line based on the vertical angle
+ d <- c(0, -1, tan(theta_deg * pi / 180))
+
+ # Calculate the endpoint of the line based on the vertical angle
+ distance_to_edge <- TopTowerSW$y - min(cloud$y)
+ z_change_due_to_angle <- distance_to_edge * d[3]
+ end_point <- c(TopTowerSW$x, TopTowerSW$y - distance_to_edge, TopTowerSW$z + z_change_due_to_angle)
+
+ return(end_point)
+}
-plot3d(Data, aspect = "iso", size = 0.2)
-pan3d(3)
-spheres3d(Devise, radius = 0.5, col = "red")
+# For demonstration: Compute endpoint with a vertical angle of 30 degrees
+theta_deg <- -60
+end_point <- compute_endpoint(TopTowerSW, cloud, theta_deg)
+segments3d(rbind(TopTowerSW, end_point), col = 'red')
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## 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")
+## spheres3d(CornA, radius = 1, col = "yellow")
## Left click on a ground point in the point cloud
## Right click to scape
diff --git a/CodeR/04_Compute_angles_cone_clear.r b/CodeR/04_Compute_angles_cone_clear.r
new file mode 100644
index 0000000000000000000000000000000000000000..38faa8ea045960c915db50211181381b99f8f6ef
--- /dev/null
+++ b/CodeR/04_Compute_angles_cone_clear.r
@@ -0,0 +1,147 @@
+rm(list=ls())
+gc()
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## packages functions
+rm(list=ls())
+
+library("rlas")
+library("rgl")
+library("stringr")
+library("VoxR")
+library("data.table")
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## DataR and CodeR
+## DataR <- paste0("/home/nibio/REPOS/holger_tower_pov/DataR")
+## DataR <- paste0("/home/pepito/Documents/holger_tower_pov/DataR")
+## DataR <- paste0("/home/nica/Documents/Holger_tower_pov/DataR")
+## CodeR <- paste0("/home/nica/Documents/Holger_tower_pov/CodeR")
+
+## DataR <- paste0("/home/nibio/REPOS/holger_tower_pov/DataR")
+## CodeR <- paste0("/home/nibio/REPOS/holger_tower_pov/CodeR")
+
+CodeR <- "/home/pepito/Documents/holger_tower_pov/CodeR/"
+DataR <- "/home/pepito/Desktop/Temp_Holger_Data/"
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## Functions
+source(paste0(CodeR, "/Functions/Functions.r"))
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## plot data - full path
+DataPaths <- list.files(DataR,
+ ## pattern = "plot_",
+ ## pattern = "Plot_",
+ pattern = "merg",
+ ## pattern = ".rds",
+ ## pattern = "Tower",
+ full.names = T)
+
+## Data <- read.las(DataPaths[1])
+Data <- fst::read_fst(DataPaths[1], as.data.table = T)
+Data <- Data[, 1:3]
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+Data <- Data[sample(1:nrow(Data), nrow(Data)*0.10), ]
+gc()
+
+Data[, X := Scale(Data$X)]
+Data[, Y := Scale(Data$Y)]
+
+## plot3d(Data, aspect = "iso", size = 0.2)
+## plot3d(Data, aspect = "iso", size = 3)
+## pan3d(3)
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## Devise location is not in the area mentioned in the email.
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## A point at the topo ofthe tower, SW edge
+TopTowerSW <- data.table(X = 173.4451,
+ Y = 104.955,
+ Z = 322.6828)
+TopTowerSE <- data.table(X = 175.3196,
+ Y = 105.2212,
+ Z = 322.7785)
+
+## spheres3d(TopTowerSE, radius = 0.3, col = "yellow")
+## spheres3d(TopTowerSW, radius = 0.3, col = "red")
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+# Function to compute the distance between points and the line
+distance2line <- function(cloud, TopTowerP, theta_deg) {
+
+ theta_rad <- theta_deg * pi / 180
+
+ d <- c(0, -1, tan(theta_rad))
+
+ cloud_south_indices <- which(cloud$Y < TopTowerP$Y)
+ cloud_south <- cloud[cloud_south_indices, ]
+
+ cloud_matrix <- as.matrix(cloud_south)
+ TopTowerP_matrix <- matrix(rep(as.matrix(TopTowerP),
+ nrow(cloud_south)),
+ ncol=3, byrow=TRUE)
+
+ r <- cloud_matrix - TopTowerP_matrix
+
+ cross_product <- cbind(r[,2]*d[3] - r[,3]*d[2],
+ r[,3]*d[1] - r[,1]*d[3],
+ r[,1]*d[2] - r[,2]*d[1])
+
+ distances <- sqrt(rowSums(cross_product^2)) / sqrt(sum(d^2))
+ Result <- data.table(distances = distances,
+ ids = cloud_south_indices)
+ gc()
+ return(Result)
+}
+
+# Example usage: Compute the distance with a vertical angle of 45 degrees
+distances <- distance2line(Data, TopTowerSW, -40)
+
+## plot3d(Data, aspect = "iso", size = 0.2)
+## spheres3d(TopTowerSW, radius = 1, col = "red")
+
+## plot3d(Data[distances[distances < 40, ]$ids, ], size = 10, col = "red", add = T)
+## pan3d(3)
+
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+Data[distances$ids, CO := distances$distances]
+
+
+## plot3d(Data, aspect = "iso", size = 0.5, col = "grey21")
+## plot3d(Data[CO < 5 , ], aspect = "iso", size = 4, col = "red", add = T)
+## spheres3d(CornA, radius = 1, col = "yellow")
+## spheres3d(CornB, radius = 1, col = "red")
+## pan3d(3)
+
+Data[!is.na(CO),
+ H := dist3d2(X, Y, Z,
+ TopTowerSW$X,
+ TopTowerSW$Y,
+ TopTowerSW$Z),
+ by = seq_len(nrow(Data))]
+
+Data[, Angle := CO/H]
+Data[, Angle := (asin(Angle))*180/pi]
+
+OpeningAngl <- 20
+
+plot3d(Data[Angle < OpeningAngl, ],
+ aspect = "iso", size = 1.5, col = "red")
+spheres3d(TopTowerSW, radius = 1, col = "yellow")
+
+plot3d(Data[Angle > OpeningAngl, ], aspect = "iso",
+ size = 0.5, col = "grey21")
+
+spheres3d(TopTowerSW, radius = 1, col = "red")
+
+plot3d(Data[Angle < OpeningAngl, ], aspect = "iso", size = 0.5, col = "red", add = T)
+pan3d(3)
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+
+
+
+
+
diff --git a/CodeR/04_Compute_angles_cone_clear.r~ b/CodeR/04_Compute_angles_cone_clear.r~
new file mode 100644
index 0000000000000000000000000000000000000000..c312c9b821d7bcd268da3617db3254c4b77d447e
--- /dev/null
+++ b/CodeR/04_Compute_angles_cone_clear.r~
@@ -0,0 +1,224 @@
+rm(list=ls())
+gc()
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## packages functions
+rm(list=ls())
+
+library("rlas")
+library("rgl")
+library("stringr")
+library("VoxR")
+library("data.table")
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## DataR and CodeR
+## DataR <- paste0("/home/nibio/REPOS/holger_tower_pov/DataR")
+## DataR <- paste0("/home/pepito/Documents/holger_tower_pov/DataR")
+## DataR <- paste0("/home/nica/Documents/Holger_tower_pov/DataR")
+## CodeR <- paste0("/home/nica/Documents/Holger_tower_pov/CodeR")
+
+## DataR <- paste0("/home/nibio/REPOS/holger_tower_pov/DataR")
+## CodeR <- paste0("/home/nibio/REPOS/holger_tower_pov/CodeR")
+
+CodeR <- "/home/pepito/Documents/holger_tower_pov/CodeR/"
+DataR <- "/home/pepito/Desktop/Temp_Holger_Data/"
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## Functions
+source(paste0(CodeR, "/Functions/Functions.r"))
+objects()
+
+dist3d
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## plot data - full path
+DataPaths <- list.files(DataR,
+ ## pattern = "plot_",
+ ## pattern = "Plot_",
+ pattern = "merg",
+ ## pattern = ".rds",
+ ## pattern = "Tower",
+ full.names = T)
+
+## Data <- read.las(DataPaths[1])
+Data <- fst::read_fst(DataPaths[1], as.data.table = T)
+Data <- Data[, 1:3]
+
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## Hi Nico,
+
+## Assuming the coordinates in the PointCloud are UTM,
+## the preferred location of the SIF sensor would be at
+
+## N 6694611.571 E 614675.0338 and height 315.8627 masl
+
+## This is a point in the middle of the SW-SE edge of the tower
+## top, and the preferred viewing direction is more or less south
+## since there is a dense group of trees there; the crucial
+## question is the viewing angle.
+
+## Cheers, Holger
+## Devise <- data.table(X = 614675.0338,
+## Y = 6694611.571,
+## Z = 315.8627)
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+Data <- Data[sample(1:nrow(Data), nrow(Data)*0.10), ]
+gc()
+
+## Data <- rbind(Data, Devise)
+
+Data[, X := Scale(Data$X)]
+Data[, Y := Scale(Data$Y)]
+## Devise <- tail(Data, 1)
+
+plot3d(Data, aspect = "iso", size = 0.2)
+pan3d(3)
+## spheres3d(Devise, radius = 0.5, col = "red")
+
+## Devise location is not in the area mentioned in the email.
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## Select the top of the tower
+## ids <- plot3d(Data, aspect = "iso", size = 4)
+## pan3d(3)
+## TopTowerIndex <- selectpoints3d(ids["data"], value = FALSE)
+
+## TopTowerIndex[, 2]
+## plot3d(Data[-TopTowerIndex[, 2], ], aspect = "iso", size = 3)
+## pan3d(3)
+## plot3d(Data[TopTowerIndex[, 2], ], aspect = "iso",
+## size = 3, col = "red", add = T)
+
+## plot3d(Data[TopTowerIndex[, 2], ], aspect = "iso",
+## size = 5, col = "red")
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## TopTowerP <- Data[TopTowerIndex[, 2], ]
+## saveRDS(TopTowerP, paste0(DataR, "TopTowerP.rds"))
+## TopTowerP <- readRDS(paste0(DataR, "TopTowerP.rds"))
+
+## plot3d(TopTowerP, aspect = "iso",
+## size = 10, col = "red")
+
+## ToWerCenter <- TopTowerP[10, ]
+## spheres3d(ToWerCenter, radius = 0.1, col = "black")
+
+## TopTowerSW <- TopTowerP[306, ]
+## spheres3d(TopTowerSW, radius = 0.1, col = "green")
+
+## TopTowerSE <- TopTowerP[168, ]
+## spheres3d(TopTowerSE, radius = 0.1, col = "blue")
+
+
+
+
+cloud <- data.frame(x = runif(10000, -20, 20), y = runif(10000, -20, 20), z = runif(10000, 0, 40))
+
+# Define the TopTowerSW point
+TopTowerSW <- data.frame(x = 0.5, y = 0.5, x= 40)
+
+# Function to compute the distance between points and the line
+distance_to_line <- function(cloud, TopTowerSW, theta_deg) {
+ # Convert angle to radians
+ theta_rad <- theta_deg * pi / 180
+
+ # Direction vector for the line
+ d <- c(0, -1, tan(theta_rad))
+
+ # Filter the cloud to retain only points to the south of TopTowerSW
+ cloud_south <- cloud[cloud$y < TopTowerSW$y, ]
+
+ # Convert both data.frames to matrices for easier computation
+ cloud_matrix <- as.matrix(cloud_south)
+ TopTowerSW_matrix <- matrix(rep(as.matrix(TopTowerSW), nrow(cloud_south)), ncol=3, byrow=TRUE)
+
+ # Compute the position vectors relative to TopTowerSW for all points
+ r <- cloud_matrix - TopTowerSW_matrix
+
+ # Compute the cross product between r and d for all points
+ cross_product <- cbind(r[,2]*d[3] - r[,3]*d[2], r[,3]*d[1] - r[,1]*d[3], r[,1]*d[2] - r[,2]*d[1])
+
+ # Compute the distance using the formula mentioned
+ distances <- sqrt(rowSums(cross_product^2)) / sqrt(sum(d^2))
+
+ return(distances)
+}
+
+# Example usage: Compute the distance with a vertical angle of 45 degrees
+distances <- distance_to_line(cloud, TopTowerSW, -60)
+
+
+plot3d(cloud, size = 5)
+plot3d(TopTowerSW, size = 10, add = T, col = "red")
+
+Puntitos <- as.numeric(names(distances[distances < 4]))
+plot3d(cloud[Puntitos, ], size = 10, col = "blue", add = T)
+
+compute_endpoint <- function(TopTowerSW, cloud, theta_deg) {
+ # Direction vector for the line based on the vertical angle
+ d <- c(0, -1, tan(theta_deg * pi / 180))
+
+ # Calculate the endpoint of the line based on the vertical angle
+ distance_to_edge <- TopTowerSW$y - min(cloud$y)
+ z_change_due_to_angle <- distance_to_edge * d[3]
+ end_point <- c(TopTowerSW$x, TopTowerSW$y - distance_to_edge, TopTowerSW$z + z_change_due_to_angle)
+
+ return(end_point)
+}
+
+# For demonstration: Compute endpoint with a vertical angle of 30 degrees
+theta_deg <- -60
+end_point <- compute_endpoint(TopTowerSW, cloud, theta_deg)
+segments3d(rbind(TopTowerSW, end_point), col = 'red')
+
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## 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()
+
+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 < 5 , ], aspect = "iso", size = 4, col = "red", add = T)
+## spheres3d(CornA, radius = 1, col = "yellow")
+## spheres3d(CornB, radius = 1, col = "red")
+## pan3d(3)
+
+Data[, H := dist3d2(X, Y, Z,
+ CornA$X, CornA$Y, CornA$Z),
+ by = seq_len(nrow(Data))]
+
+Data[, Angle := CO/H]
+Data[, Angle := (asin(Angle))*180/pi]
+
+OpeningAngl <- 20
+
+plot3d(Data[Angle < OpeningAngl, ], aspect = "iso", size = 1.5, col = "red")
+spheres3d(CornA, radius = 1, col = "yellow")
+spheres3d(CornB, radius = 1, col = "red")
+
+plot3d(Data[Angle > OpeningAngl, ], aspect = "iso", size = 0.5, col = "grey21")
+spheres3d(CornA, radius = 1, col = "yellow")
+spheres3d(CornB, radius = 1, col = "red")
+plot3d(Data[Angle < OpeningAngl, ], aspect = "iso", size = 0.5, col = "red", add = T)
+pan3d(3)
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+
+
+
+
+