<p>多亏了卡蒂娅,她完成了这项工作!</strong></p>
<p>要计算三维空间中给定点列表的半径,我使用以下代码:</p>
<pre><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)
</code></pre>
<p>再次感谢你!代码很好!在</p>