From 37ab0cb089057d8875897b8e7ca8d99dcec641b1 Mon Sep 17 00:00:00 2001
From: nicsoattc <nicolas.cattaneo@nibio.no>
Date: Mon, 23 Oct 2023 14:04:42 +0200
Subject: [PATCH] update repo
---
CodeR/04_Compute_angles_cone_clear.r | 34 ++--
CodeR/05_Compute_angles_cone_MAKEA_FUNCTION.r | 143 +++++++++++++++++
.../05_Compute_angles_cone_MAKEA_FUNCTION.r~ | 149 ++++++++++++++++++
MASTER.zip | Bin 0 -> 3598 bytes
MASTER/00_apply_sifFootP.r | 68 ++++++++
MASTER/00_apply_sifFootP.r~ | 68 ++++++++
MASTER/FUNCTIONS/pan3d.r | 26 +++
MASTER/FUNCTIONS/sifFootP.r | 49 ++++++
MASTER/FUNCTIONS/sifFootP.r~ | 45 ++++++
9 files changed, 566 insertions(+), 16 deletions(-)
create mode 100644 CodeR/05_Compute_angles_cone_MAKEA_FUNCTION.r
create mode 100644 CodeR/05_Compute_angles_cone_MAKEA_FUNCTION.r~
create mode 100644 MASTER.zip
create mode 100644 MASTER/00_apply_sifFootP.r
create mode 100644 MASTER/00_apply_sifFootP.r~
create mode 100644 MASTER/FUNCTIONS/pan3d.r
create mode 100644 MASTER/FUNCTIONS/sifFootP.r
create mode 100644 MASTER/FUNCTIONS/sifFootP.r~
diff --git a/CodeR/04_Compute_angles_cone_clear.r b/CodeR/04_Compute_angles_cone_clear.r
index 38faa8e..7e7693d 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 0000000..9a8e198
--- /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 0000000..7e7693d
--- /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
GIT binary patch
literal 3598
zcmWIWW@Zs#0D+xB;o%?}hS?bu7<?UrLtKOOLqm8O*yVN2Qd#xPQbBZS1vdjD%L`@(
z1~36JgoA;D;djowAeaFlS{h-9fkAv?K|xMsd~s%)TYi2?fL>7@*x>Fr-3fbD%~HkH
z%u>+|W@ZTRX6FDI!~`N37#I*%A-i`kI|G9hUaP=1ft(Jq1w^CSv^F$1|CWKkzj^V6
zUf&g0W;87M!pfT|`gTigMmN{elHkytCj}Q-7N%G)&AC;7-;&X3`$qAt{}{Htm%D%W
z@pGA(RyKjVj!%8OD@?{}UFzF62X>3|7CW6gq_5+~dRd9PRLLmyWco&(c{Uzu;(fZ7
z!P|Ck?ECpFD061Y#yfYUx|gp$TYALgvq|-XjQw8!o_Fu-_`r4VUE|$3og%X%g*H7~
z$<Zxhympp6OMl?KY{TXEV{FXM`?oFLy!zSi*WLeKJbrrg!ZKq&&DH-c7k}9v-&wC0
zvE|@ip;n%YCud%|Q8ewVsmA6xsg6$scfTu>pXzg;dm;N&>5#`$Qx3JV71wTXzggtB
z>!R^#1H%QKMsFm~R$OjbTr_QtS(xP?BaXF!X&0Dpe}C^RB=&Tkf~bUMwe=a5O}v>$
zcFJZJeKVRi=|Z)E(#M%BzSn}6?vCx$*|mJ;X1x`RK9yTzv;|t@5AqypStI!Fm&fP-
zI~m{R6@+&^iBdEO&^h~3LD%7FhoFdMi_%%{v#Vy;x2>Bqx6FKUQrxiu>45X)(Q&(7
zt;@f98@TOkh!i^iVB*K7QzF4L7)~qgXq7OY;-mZ|D8O)K(@!4bRf^pgzN|U;TGC0t
zQ76*Q+v%y}{3HFF?wog>l*zGELn&LrXw@Ns?N`<q9}D1qS*N|Te46(p&sl47^EU62
zP*3$=y?3wg?%lih?Obbjf3>%An)-#ejx#+xkI(L`-fndFU&=1k;FIT1ZsJ;YdeW8?
zyv|prNbP1kEPX)3ZG)SV{7l_j3yiXPErf#{e;a)G_TBpJTd&0;C)RF#FI#R^|KXQQ
zT;uy{cK3f1qmSGC&OOs7zNMf4pVa2ZKkIJ)Z(q>f?mpMcY}xwe6;Doi*r~7c%~tga
zFrHXz8S?4s7l#W9XNw*Sc*p$6-Tj7h>iK|0i;`a%Nhw555m|OSDvmc=;d#<N(X!2~
zYKE%}o@&2f(f?Z^QL;bBuI1p<8&S4%&K<~DlHcN&=pB@06duWUCOJSn%6h^K6YHZh
z-c;9}-g@oiKXCT)=6Sa9I}-x~C=)@l7bveiLn{sh5!uTv)XzD@)88)`QZN+inx(Sn
znWchgv|I%-4pcB88~2fwfdQ#tKp2N88Nk^UWMH^1R)bqZ?t`*x-ShPyH2BT9?S)J-
zXFW(<w0f6%!09DHQx*yAPT+iVrtC`D<$vGj&J~c!U{q|Wy<h(Q?$387_4_7X^0~O7
zd0z2EDLvMCUV3LLw2P${%si9!sdZ&$!<n|3PiCz>d7sVcgX__=L8@0SeNVgMCh2Y~
zQuo$+g7LB)@(cVp{Dc~{r4;!d)&)uDuT0G|^VZYwpSkVGbY6e+d5f;e@bWK;DlYQb
zwuHInW8ZCtqsJ4=I|Kh*4HCYq%XcGL%yHEUxjjd;raffd%rxnMpngruB>k-!(=(N(
z-oB{du-g1VUt0P8x{9eG>mJm_%@uwd$2jFm)l5y#&sQv-alLA}s?_nwES4uJ;c#LC
ze_&16O*ye?b2bzdbP9{D<Ms=&_;6`qjPvCwZx_^SH!D8>U9-H<_gAnBN6B2~7d7oG
z_ga)X2v08Y-E_;SW{142!E^8PD_Qqn4qkGxz-9est?9Z6tQj_r_xD(=eHCAG=+BYT
zaFq%h%iNae*GHaq$p1;twVkk>qxN5f&AjUmK6`&SF7-$7!I9Zg{y7c%y1!0~jD76+
z-pnq!UGO2(Ze_*X8~l5Frr*5(f6vSmpMs6s(tZbAT(BV}I+IN#{_t*AgSnE;da?`E
zu5XH5-S@0sLgLq5ea5Dj8)j@x$`8Du(yBZ+cSe#_O{F#OG=cM9HTyf&J|3>-6n)9l
zuW~qO<w?t+UDL&@5=9n&6%q1!SQqwc?=GkE9zS(~2zlMD{Hkv{yJp64#(X?onDr<7
z4?CzT(VMzS(1eMB0h9z0i8I>|QT4rHW?+!Omevar^NdsUiohxJx2jp{H8rzTPRuH6
zZSd*5W&?pepGCE=c&IgP;e4xW`E`!IN3p>Jm%Lc%n`#ZK{#5_v=ojG(&8ppKI`8J)
zkL$e8TxQzFIY(<=wwA>yt+lNE540IJOPp1irpxtoX6wY}DaV4MITw2=_oqJb`58My
zt!*3krPXp;u21B$J*96*m35{@Ehuov*JL}`I%$`P*?oui8#p7M&X{O;nNc`O^Qd)U
zsMWf{ZBhy*iQHE&XPx2>%+J}&G3)qmA^w+!@3wx<=H=XC|MdKEiFslA`+E(GHqMs#
zthTQ4_Nq+>W}PWpo_J%c$luAQ?WN{?XFF-)bbob*<*e2P>UCVn%V+Ol*>WUC^1Rv{
z_KCHPhi`?v;G6w!f?Rv$2WL~Mhj&X>iJxZ8FgYc5FjoJ~Z>>e&r+w&>JYDc_VP3u7
zxBD9VOwYV7S`&D~_+ZAtJze50h9##0=H%pPoM?S^B<6-orD|ZxyL&r+GGBin4T`(n
z0h&+47#SEqF@%V_hyX;~&17X@kiiyrh=K*u<XEX;mby*DEEQC`pcPrIA-D6S3<YYR
zYsVaE_C1>_?z&2=Vui=zFPv8dG+A;lY$`h-bR+t<_N<V<-~0B)?3uvi#8G|DvN+9V
zo8bC#ubYOiwqM9rG}fMSBw2Nd?7Gl{8naf^EZP>+uenJ(?{f4$mHI4>gW^-Ii{~1-
z-)|2Te#%rD@MY<1N4psh^<39-xCC|bn>g7V`k$hA|IUg->z+9oO)PxeaeL~Cu*Sv5
zM85VflL{{KDqEsl^Q}?G^<?HP->j`a-mX}gIh|X^;IvBn8pfLTt2W2opD}u@P>=iR
zv3A?koa>oJ!LJ?r7af%@V4pm{|9-Py*8L9~)?J(VeF0;Xr1AN!Nv>(O#f-Iws@xtg
zIhnSX<;I2=J6-3*lx%Bvwe}0)ih1m^wB%%`cjE#XH6EU`c`|48a}BTj3CQ7D7g68I
z*|_pELvW5kvGj8$ovoR-mvx&SKmPZiVDZkcO_m?nja{=${dv^hE8Judb3H8hOX^9c
z%s<uiGyaFRq;|D0{_MD;Txp-%9}&6i$DMnQy!ib8?Bh&UyS;4n>kiCk^lOS1p0V`B
z?3i^sKby*wzMITlxI@V)X}`)6&iBsh;=5nT*WKwY(rb;`wZVGDj&Qe6TQ}LZ<v+Zm
zap&Be!&CSIVkJV?mmKdX=WERS-KMldmT&WwKe~OdmWUYr?8?q_(5|@rKv1aF@_AvT
zdg8AI3pUO9*O^)oqE+TM=iOUZN2#-sTU(|W-k9HTfB)8cP!hbuZ6_tk#J~VbXaU}g
zOd`y<+s>d$6omgef=HrU(hx&IRVJ)04Qu*>s!WhkAiSl~m4N})JSV4R4RbxnB2Y^j
zxe)+rNrNl_;Vq3D@LB{e4ngf{9PMhDJs`7?+W;WGAcF#k0^`4q{48(}AzB1r<3a5J
zoDB>_c!7*ZF5gkZD-N6Sh;jjL5vWK2S%_9F!2AIUuWY~YaOA=WR9%2Q!N9<<rSUf#
x7Jm>~c_7;YD%g<I8>nDI*s`4ii!BK65G>{byjj^mrm-`yGwfqxV6fr_@c{lYP(J_w
literal 0
HcmV?d00001
diff --git a/MASTER/00_apply_sifFootP.r b/MASTER/00_apply_sifFootP.r
new file mode 100644
index 0000000..7998f21
--- /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 0000000..21f49dc
--- /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 0000000..03e46ca
--- /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 0000000..42fcc47
--- /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 0000000..2473905
--- /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)
+
+}
--
GitLab