diff --git a/CodeR/04_Compute_angles_cone_clear.r b/CodeR/04_Compute_angles_cone_clear.r
index 38faa8ea045960c915db50211181381b99f8f6ef..7e7693d1bb206536b4d3bed0e07abcc7bf2d791f 100644
--- a/CodeR/04_Compute_angles_cone_clear.r
+++ b/CodeR/04_Compute_angles_cone_clear.r
@@ -21,8 +21,8 @@ library("data.table")
## 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/"
+CodeR <- "/home/nica/Documents/Holger_tower_pov/CodeR/"
+DataR <- "/home/nica/Documents/Holger_tower_pov/DataR"
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## Functions
@@ -39,7 +39,7 @@ DataPaths <- list.files(DataR,
full.names = T)
## Data <- read.las(DataPaths[1])
-Data <- fst::read_fst(DataPaths[1], as.data.table = T)
+Data <- fst::read_fst(DataPaths[2], as.data.table = T)
Data <- Data[, 1:3]
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@@ -109,33 +109,35 @@ distances <- distance2line(Data, TopTowerSW, -40)
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
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[!is.na(CO),
- H := dist3d2(X, Y, Z,
- TopTowerSW$X,
- TopTowerSW$Y,
- TopTowerSW$Z),
- by = seq_len(nrow(Data))]
+ H := sqrt((X - TopTowerSW$X)^2 +
+ (Y - TopTowerSW$Y)^2 +
+ (Z - TopTowerSW$Z)^2)]
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",
+plot3d(Data, aspect = "iso",
size = 0.5, col = "grey21")
-
-spheres3d(TopTowerSW, radius = 1, col = "red")
+spheres3d(TopTowerSW, radius = 1, col = "yellow")
+plot3d(Data[Angle < OpeningAngl, ],
+ aspect = "iso", size = 1.5, col = "red",
+ add = T)
plot3d(Data[Angle < OpeningAngl, ], aspect = "iso", size = 0.5, col = "red", add = T)
pan3d(3)
diff --git a/CodeR/05_Compute_angles_cone_MAKEA_FUNCTION.r b/CodeR/05_Compute_angles_cone_MAKEA_FUNCTION.r
new file mode 100644
index 0000000000000000000000000000000000000000..9a8e198f563c59e71ba16f597cc7a590915606f3
--- /dev/null
+++ b/CodeR/05_Compute_angles_cone_MAKEA_FUNCTION.r
@@ -0,0 +1,143 @@
+rm(list=ls())
+gc()
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## packages functions
+rm(list=ls())
+
+library("rlas")
+library("rgl")
+library("stringr")
+library("VoxR")
+library("data.table")
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+
+CodeR <- "/home/nica/Documents/Holger_tower_pov/CodeR/"
+DataR <- "/home/nica/Documents/Holger_tower_pov/DataR"
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## pan3d function
+pan3d <- function(button) {
+ library("rgl")
+ start <- list()
+ begin <- function(x, y) {
+ start$userMatrix <<- par3d("userMatrix")
+ start$viewport <<- par3d("viewport")
+ start$scale <<- par3d("scale")
+ start$projection <<- rgl.projection()
+ start$pos <<- rgl.window2user( x/start$viewport[3],
+ 1 - y/start$viewport[4], 0.5,
+ projection = start$projection)
+ }
+
+ update <- function(x, y) {
+ xlat <- (rgl.window2user( x/start$viewport[3],
+ 1 - y/start$viewport[4], 0.5,
+ projection = start$projection) -
+ start$pos)*start$scale
+
+ mouseMatrix <- translationMatrix(xlat[1], xlat[2], xlat[3])
+ par3d(userMatrix = start$userMatrix %*% t(mouseMatrix) )
+ }
+
+ rgl.setMouseCallbacks(button, begin, update)
+ cat("Callbacks set on button", button, "of rgl device", rgl.cur(), "\n")
+}
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## plot data - full path
+DataPaths <- list.files(DataR,
+ pattern = "22_classified.rds",
+ ## pattern = "22_classified.las",
+ full.names = T)
+
+Data <- fst::read_fst(DataPaths, as.data.table = T)
+## Data <- read.las(DataPaths)
+Data <- Data[, 1:3]
+
+## Use a sample to speed up computations
+Data <- Data[sample(1:nrow(Data), nrow(Data)*0.10), ]
+gc()
+
+## Very high values of X and Y generate problems when rendering
+## the 3d graphics. Center the point cloud
+Scale <- function(Vector){
+ ((Vector - min(Vector))/
+ (max(Vector)-min(Vector)))*
+ (max(Vector)-min(Vector))}
+
+Data[, X := Scale(Data$X)]
+Data[, Y := Scale(Data$Y)]
+
+
+## A point at the top of the 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)
+
+## plot3d(Data, aspect = "iso",
+## size = 3, col = "grey21")
+## spheres3d(TopTowerSE, radius = 0.3, col = "yellow")
+## spheres3d(TopTowerSW, radius = 0.3, col = "red")
+## pan3d(3)
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+# Function to compute the distance between points and the line
+
+sifFootP <- function(PointCloud, TopTowerP, thetaDeg){
+
+ theta_rad <- thetaDeg * pi / 180
+ d <- c(0, -1, tan(theta_rad))
+
+ cloud_south_indices <- which(PointCloud$Y < TopTowerP$Y)
+ cloud_south <- PointCloud[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))
+ distances <- data.table(distances = distances,
+ ids = cloud_south_indices)
+
+
+ PointCloud[distances$ids, CO := distances$distances]
+
+ PointCloud[!is.na(CO),
+ H := sqrt((X - TopTowerP$X)^2 +
+ (Y - TopTowerP$Y)^2 +
+ (Z - TopTowerP$Z)^2)]
+
+ PointCloud[, Angle := CO/H]
+ PointCloud[, Angle := (asin(Angle))*180/pi]
+
+ SIFOpeningAngl <- 20
+
+ plot3d(Data, aspect = "iso",
+ size = 0.5, col = "grey21")
+
+ spheres3d(TopTowerP, radius = 0.5, col = "green")
+ plot3d(PointCloud[Angle < SIFOpeningAngl, ],
+ aspect = "iso", size = 1.5, col = "red",
+ add = T)
+
+}
+
+sifFootP(PointCloud = Data,
+ TopTowerP = TopTowerSW,
+ thetaDeg = -30)
+
+
+
+
diff --git a/CodeR/05_Compute_angles_cone_MAKEA_FUNCTION.r~ b/CodeR/05_Compute_angles_cone_MAKEA_FUNCTION.r~
new file mode 100644
index 0000000000000000000000000000000000000000..7e7693d1bb206536b4d3bed0e07abcc7bf2d791f
--- /dev/null
+++ b/CodeR/05_Compute_angles_cone_MAKEA_FUNCTION.r~
@@ -0,0 +1,149 @@
+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/nica/Documents/Holger_tower_pov/CodeR/"
+DataR <- "/home/nica/Documents/Holger_tower_pov/DataR"
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## 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[2], 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[!is.na(CO),
+ H := sqrt((X - TopTowerSW$X)^2 +
+ (Y - TopTowerSW$Y)^2 +
+ (Z - TopTowerSW$Z)^2)]
+
+Data[, Angle := CO/H]
+Data[, Angle := (asin(Angle))*180/pi]
+
+OpeningAngl <- 20
+
+plot3d(Data, aspect = "iso",
+ size = 0.5, col = "grey21")
+spheres3d(TopTowerSW, radius = 1, col = "yellow")
+plot3d(Data[Angle < OpeningAngl, ],
+ aspect = "iso", size = 1.5, col = "red",
+ add = T)
+
+plot3d(Data[Angle < OpeningAngl, ], aspect = "iso", size = 0.5, col = "red", add = T)
+pan3d(3)
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+
+
+
+
+
diff --git a/MASTER.zip b/MASTER.zip
new file mode 100644
index 0000000000000000000000000000000000000000..1df67f9d9b55b35afaa4cad081955d0e5211bfc4
Binary files /dev/null and b/MASTER.zip differ
diff --git a/MASTER/00_apply_sifFootP.r b/MASTER/00_apply_sifFootP.r
new file mode 100644
index 0000000000000000000000000000000000000000..7998f2121fb5d406127e1384d44ddd9b14df9b9a
--- /dev/null
+++ b/MASTER/00_apply_sifFootP.r
@@ -0,0 +1,68 @@
+rm(list=ls())
+gc()
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## packages functions
+rm(list=ls())
+
+library("rlas")
+library("rgl")
+library("stringr")
+library("VoxR")
+library("data.table")
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+MyCodeR <- "/home/nica/Documents/Holger_tower_pov/MASTER/FUNCTIONS"
+MyDataR <- "/home/nica/Documents/Holger_tower_pov/DataR"
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## load functions
+source(paste0(MyCodeR, "/pan3d.r"))
+source(paste0(MyCodeR, "/sifFootP.r"))
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## plot data - full path
+DataPaths <- list.files(MyDataR,
+ pattern = "22_classified.rds",
+ ## pattern = "22_classified.las",
+ full.names = T)
+
+Data <- fst::read_fst(DataPaths, as.data.table = T)
+## Data <- read.las(DataPaths)
+Data <- Data[, 1:3]
+
+## Use a sample to speed up computations
+Data <- Data[sample(1:nrow(Data), nrow(Data)*0.10), ]
+gc()
+
+## Very high values of X and Y generate problems when rendering
+## the 3d graphics. Center the point cloud
+Scale <- function(Vector){
+ ((Vector - min(Vector))/
+ (max(Vector)-min(Vector)))*
+ (max(Vector)-min(Vector))}
+
+Data[, X := Scale(Data$X)]
+Data[, Y := Scale(Data$Y)]
+
+
+## A point at the top of the 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)
+
+## plot3d(Data, aspect = "iso",
+## size = 3, col = "grey21")
+## spheres3d(TopTowerSE, radius = 0.3, col = "yellow")
+## spheres3d(TopTowerSW, radius = 0.3, col = "red")
+## pan3d(3)
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## Compute POV
+sifFootP(PointCloud = Data,
+ TopTowerP = TopTowerSE,
+ thetaDeg = -45)
+pan3d(3)
diff --git a/MASTER/00_apply_sifFootP.r~ b/MASTER/00_apply_sifFootP.r~
new file mode 100644
index 0000000000000000000000000000000000000000..21f49dc60993496097cf53d1eaaad116e1649227
--- /dev/null
+++ b/MASTER/00_apply_sifFootP.r~
@@ -0,0 +1,68 @@
+rm(list=ls())
+gc()
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## packages functions
+rm(list=ls())
+
+library("rlas")
+library("rgl")
+library("stringr")
+library("VoxR")
+library("data.table")
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+MyCodeR <- "/home/nica/Documents/Holger_tower_pov/MASTER/FUNCTIONS"
+MyDataR <- "/home/nica/Documents/Holger_tower_pov/DataR"
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## load functions
+source(paste0(MyCodeR, "/pan3d.r"))
+source(paste0(MyCodeR, "/sifFootP.r"))
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## plot data - full path
+DataPaths <- list.files(MyDataR,
+ pattern = "22_classified.rds",
+ ## pattern = "22_classified.las",
+ full.names = T)
+
+Data <- fst::read_fst(DataPaths, as.data.table = T)
+## Data <- read.las(DataPaths)
+Data <- Data[, 1:3]
+
+## Use a sample to speed up computations
+Data <- Data[sample(1:nrow(Data), nrow(Data)*0.10), ]
+gc()
+
+## Very high values of X and Y generate problems when rendering
+## the 3d graphics. Center the point cloud
+Scale <- function(Vector){
+ ((Vector - min(Vector))/
+ (max(Vector)-min(Vector)))*
+ (max(Vector)-min(Vector))}
+
+Data[, X := Scale(Data$X)]
+Data[, Y := Scale(Data$Y)]
+
+
+## A point at the top of the 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)
+
+## plot3d(Data, aspect = "iso",
+## size = 3, col = "grey21")
+## spheres3d(TopTowerSE, radius = 0.3, col = "yellow")
+## spheres3d(TopTowerSW, radius = 0.3, col = "red")
+## pan3d(3)
+
+## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+## Compute POV
+sifFootP(PointCloud = Data,
+ TopTowerP = TopTowerSE,
+ thetaDeg = -20)
+pan3d(3)
diff --git a/MASTER/FUNCTIONS/pan3d.r b/MASTER/FUNCTIONS/pan3d.r
new file mode 100644
index 0000000000000000000000000000000000000000..03e46ca3f18581e5fafcc8f3e9a21b675cd35dbc
--- /dev/null
+++ b/MASTER/FUNCTIONS/pan3d.r
@@ -0,0 +1,26 @@
+pan3d <- function(button) {
+ library("rgl")
+ start <- list()
+ begin <- function(x, y) {
+ start$userMatrix <<- par3d("userMatrix")
+ start$viewport <<- par3d("viewport")
+ start$scale <<- par3d("scale")
+ start$projection <<- rgl.projection()
+ start$pos <<- rgl.window2user( x/start$viewport[3],
+ 1 - y/start$viewport[4], 0.5,
+ projection = start$projection)
+ }
+
+ update <- function(x, y) {
+ xlat <- (rgl.window2user( x/start$viewport[3],
+ 1 - y/start$viewport[4], 0.5,
+ projection = start$projection) -
+ start$pos)*start$scale
+
+ mouseMatrix <- translationMatrix(xlat[1], xlat[2], xlat[3])
+ par3d(userMatrix = start$userMatrix %*% t(mouseMatrix) )
+ }
+
+ rgl.setMouseCallbacks(button, begin, update)
+ cat("Callbacks set on button", button, "of rgl device", rgl.cur(), "\n")
+}
diff --git a/MASTER/FUNCTIONS/sifFootP.r b/MASTER/FUNCTIONS/sifFootP.r
new file mode 100644
index 0000000000000000000000000000000000000000..42fcc47a96e8ce12e5328948e3bed230631ed08b
--- /dev/null
+++ b/MASTER/FUNCTIONS/sifFootP.r
@@ -0,0 +1,49 @@
+sifFootP <- function(PointCloud, TopTowerP, thetaDeg){
+
+ theta_rad <- thetaDeg * pi / 180
+ d <- c(0, -1, tan(theta_rad))
+
+ cloud_south_indices <- which(PointCloud$Y < TopTowerP$Y)
+ cloud_south <- PointCloud[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))
+ distances <- data.table(distances = distances,
+ ids = cloud_south_indices)
+
+
+ PointCloud[distances$ids, CO := distances$distances]
+
+ PointCloud[!is.na(CO),
+ H := sqrt((X - TopTowerP$X)^2 +
+ (Y - TopTowerP$Y)^2 +
+ (Z - TopTowerP$Z)^2)]
+
+ PointCloud[, Angle := CO/H]
+ PointCloud[, Angle := (asin(Angle))*180/pi]
+
+ SIFOpeningAngl <- 23
+
+ plot3d(Data, aspect = "iso",
+ size = 0.5, col = "grey21")
+
+ spheres3d(TopTowerP, radius = 0.5, col = "green")
+ plot3d(PointCloud[Angle < SIFOpeningAngl, ],
+ aspect = "iso", size = 1.5, col = "red",
+ add = T)
+
+ ## RedPOints <- PointCloud[Angle < SIFOpeningAngl, ]
+ ## return(RedPOints)
+
+
+}
diff --git a/MASTER/FUNCTIONS/sifFootP.r~ b/MASTER/FUNCTIONS/sifFootP.r~
new file mode 100644
index 0000000000000000000000000000000000000000..24739058fb4a5d998a1ebeb631e53a1b01bcd295
--- /dev/null
+++ b/MASTER/FUNCTIONS/sifFootP.r~
@@ -0,0 +1,45 @@
+sifFootP <- function(PointCloud, TopTowerP, thetaDeg){
+
+ theta_rad <- thetaDeg * pi / 180
+ d <- c(0, -1, tan(theta_rad))
+
+ cloud_south_indices <- which(PointCloud$Y < TopTowerP$Y)
+ cloud_south <- PointCloud[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))
+ distances <- data.table(distances = distances,
+ ids = cloud_south_indices)
+
+
+ PointCloud[distances$ids, CO := distances$distances]
+
+ PointCloud[!is.na(CO),
+ H := sqrt((X - TopTowerP$X)^2 +
+ (Y - TopTowerP$Y)^2 +
+ (Z - TopTowerP$Z)^2)]
+
+ PointCloud[, Angle := CO/H]
+ PointCloud[, Angle := (asin(Angle))*180/pi]
+
+ SIFOpeningAngl <- 23
+
+ plot3d(Data, aspect = "iso",
+ size = 0.5, col = "grey21")
+
+ spheres3d(TopTowerP, radius = 0.5, col = "green")
+ plot3d(PointCloud[Angle < SIFOpeningAngl, ],
+ aspect = "iso", size = 1.5, col = "red",
+ add = T)
+
+}