04_Compute_angles_cone_clear.r

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)
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@