Radius of three points out of a points list for python or R

Question

For the purpose of a mathematical model of riding along an arbitrary path, I need to know the radius of each segment. My path is defined by points with the distance /dx in 3d space and listed in a text file (x;y;z).

To calculate the radius I need at least three points. Therefore, I'm looking for a script (python or R) that takes the first three points, calculate the radius and stores the value to the middle point. For the next iteration it should take the second point as start and the following two points -> r(P1,P2,P3) -> P1(x,y,z) = n; P2(x,y,z) = n + 1; P3(x,y,z) = n + 2;

Thanks for any advise and your support!

What do you mean by radius? The radius of the smallest circle containing all three points or something? For example, what if your three points are (0,0,0), (1,0,0), (2,0,0). What is the radius there? — n3utrino, May 16 '18 at 17:42
consider an arbitrary path with four points e.g. (0,0,0), (1,0,0), (2,0,0) & (2,1,0). for the first three, the radius is infinity and for the next tripple it about 0.7 — google Nutzer, May 16 '18 at 17:53
You should be able to get the code from this answer: https://stackoverflow.com/a/41145486/3080953 — c2huc2hu, May 16 '18 at 18:04

score 0 · Answer 1 · answered May 17 '18 at 03:24

# distance between 2 points
eucl.dist <- function(x1, x2) sqrt(sum((x1 - x2) ^ 2))

# area of triangle
area <- function(x1,x2,x3) {
  v1 = x1 - x2
  v2 = x1 - x3
  0.5* sqrt( (v1[2]*v2[3] - v1[3]*v2[2])**2 + 
             (v1[3]*v2[1] - v1[1]*v2[3])**2 + 
             (v1[1]*v2[2] - v1[2]*v2[1])**2)
}

# function to calculate curvature (optional for this problem)
curvature <- function(x1, x2, x3){
  4 * area(x1, x2, x3) / (eucl.dist(x1,x2)*eucl.dist(x2,x3)*eucl.dist(x3,x1)) 
}

# function to calculate radius
radius <- function(x1, x2, x3){
  0.25 *  (eucl.dist(x1,x2)*eucl.dist(x2,x3)*eucl.dist(x3,x1))  / area(x1, x2, x3) 
}


# Test on some simple 2D case  
path2D <- cbind( c(0, 1, 2, 2, 3, 4, 4, 4),
                    c(0, 1, 0,-1,-1, 0, 1, 3),
                    rep(0,8) )

names(path2D) <- x("x","y","z")

#plot the path
plot( path2D[,1], path2D[,2], type="b", 
      lwd=2, col="slateblue", 
      xlim=c(-1,5), ylim=c(-2,4),
      xlab="x", ylab="y")

# calculate vector with radii
rads <- sapply( 1 : (nrow(path2D) - 2), FUN=function(i){ radius(x1 = path2D[i,], x2 = path2D[i+1,], x3 = path2D[i+2,] )} )
# [1] 1.0000000 1.5811388 0.7071068 1.5811388 1.5811388     Inf

# Add text to the plot
text(path2D[2,1],path2D[2,2], round(rads[1],2), pos=3, cex=.8)
text(path2D[3,1],path2D[3,2], round(rads[2],2), pos=4, cex=.8)
text(path2D[4,1],path2D[4,2], round(rads[3],2), pos=2, cex=.8)
text(path2D[5,1],path2D[5,2], round(rads[4],2), pos=1, cex=.8)
text(path2D[6,1],path2D[6,2], round(rads[5],2), pos=4, cex=.8)
text(path2D[7,1],path2D[7,2], round(rads[6],2), pos=4, cex=.8)

@googleNutzer Would you like me to add more comments to the code so it would be easier to understand what happens in each step? — Katia, May 17 '18 at 15:35

score 0 · Answer 2 · answered May 17 '18 at 18:25

Many thanks to Katia, she've done the work!

To calculate the radius for a given points list in 3d space, I use the following code:

rm(list=ls()) 

setwd("...")

## ask Katia @ stackoverflow

# distance between 2 points
eucl.dist <- function(x1, x2) sqrt(sum((x1 - x2) ^ 2))

# area of triangle
area <- function(x1,x2,x3) {
  v1 = x1 - x2
  v2 = x1 - x3
  0.5* sqrt( (v1[2]*v2[3] - v1[3]*v2[2])**2 + 
               (v1[3]*v2[1] - v1[1]*v2[3])**2 + 
               (v1[1]*v2[2] - v1[2]*v2[1])**2)
}

# function to calculate curvature (optional for this problem)
curvature <- function(x1, x2, x3){
  4 * area(x1, x2, x3) / (eucl.dist(x1,x2)*eucl.dist(x2,x3)*eucl.dist(x3,x1)) 
}

# function to calculate radius
radius <- function(x1, x2, x3){
  0.25 *  (eucl.dist(x1,x2)*eucl.dist(x2,x3)*eucl.dist(x3,x1))  / area(x1, x2, x3) 
}


# reading data 
path3D <- read.table("...txt", dec=".", sep=";", header=TRUE)

colnames(path3D) <- c("x","y","z")


# calculate vector with radii
rads <- sapply( 1 : (nrow(path3D) - 2), FUN=function(i){ radius(x1 = path3D[i,], x2 = path3D[i+1,], x3 = path3D[i+2,] )} )
rads <- data.frame(rads)
#add header
colnames(rads) <- c("radius")

#considering, that for the first value no radius could be calculated
rads <- rbind( Inf, rads)
#...also for the last one
rads <- rbind(rads, Inf)
#merge the data
path3D <- cbind(path3D, data.frame(rads))

#write to file
write.table(path3D, "....txt", sep=";",row.names=FALSE)

Thank you again! Well coded!

Radius of three points out of a points list for python or R

2 Answers2