Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Torus passing by three points in rgl
# see https://laustep.github.io/stlahblog/posts/rglTorus.html
#### transformation matrix and radius ####
# plane passing by points p1, p2, p3 #
plane3pts <- function(p1,p2,p3){
xcoef = (p1[2]-p2[2])*(p2[3]-p3[3])-(p1[3]-p2[3])*(p2[2]-p3[2]);
ycoef = (p1[3]-p2[3])*(p2[1]-p3[1])-(p1[1]-p2[1])*(p2[3]-p3[3]);
zcoef = (p1[1]-p2[1])*(p2[2]-p3[2])-(p1[2]-p2[2])*(p2[1]-p3[1]);
offset = p1[1]*xcoef + p1[2]*ycoef + p1[3]*zcoef;
c(xcoef, ycoef, zcoef, offset)
}
# center, radius and normal of the circle passing by three points #
circleCenterAndRadius <- function(p1,p2,p3){
p12 <- (p1+p2)/2
p23 <- (p2+p3)/2
v12 <- p2-p1
v23 <- p3-p2
plane1 <- plane3pts(p1,p2,p3)
A <- rbind(plane1[1:3],v12,v23)
b <- c(plane1[4], sum(p12*v12), sum(p23*v23))
center <- c(solve(A) %*% b)
r <- sqrt(c(crossprod(p1-center)))
list(center=center, radius=r, normal=plane1[1:3])
}
# transformation matrix #
transfoMatrix <- function(p1,p2,p3){
crn <- circleCenterAndRadius(p1,p2,p3)
center <- crn$center
radius <- crn$radius
normal <- crn$normal
if(normal[1]==0 && normal[2]==0){
m <- rbind(cbind(diag(3), center), c(0,0,0,1))
}else{
measure <- sqrt(normal[1]^2+normal[2]^2+normal[3]^2)
normal <- normal/measure
s <- sqrt(normal[1]^2+normal[2]^2)
u <- c(normal[2]/s, -normal[1]/s, 0)
v <- geometry::extprod3d(normal, u)
m <- rbind(cbind(u, v, normal, center), c(0,0,0,1))
}
list(matrix=t(m), radius=radius)
}
#### torus as tmesh3d ####
# R: major radius
# r: minor radius
# S: number of segments for the major ring
# s: number of segments for the minor ring
# arc: the arc to draw
library(rgl)
createTorusMesh <- function(R=1, r=0.25, S=64, s=32, arc=2*pi){
vertices <- indices <- NULL
for (j in 0:s) {
v <- j / s * 2*pi
for (i in 0:S) {
u <- i / S * arc
vertex <- c(
(R + r*cos(v)) * cos(u),
(R + r*cos(v)) * sin(u),
r * sin(v)
)
vertices <- cbind(vertices, vertex)
}
}
for (j in 1:s) {
for (i in 1:S) {
a <- (S + 1) * j + i
b <- (S + 1) * (j - 1) + i
indices <- cbind(indices, c(a,b,a+1), c(b,b+1,a+1))
}
}
tmesh3d(rbind(vertices,1), indices)
}
# #### torus passing by three points ####
# library(rgl)
#
# p1 = c(1,-1,1)
# p2 = c(2,1,2)
# p3 = c(2,-2,-3)
# mr <- transfoMatrix(p1,p2,p3)
#
# tmesh0 <- createTorusMesh(R=mr$radius, r=0.1)
#
# tmesh <- transform3d(tmesh0, mr$matrix)
# shade3d(tmesh, color="blue")
#
# triangles3d(rbind(p1,p2,p3), color="red")
torusMesh3points <- function(p1, p2, p3, r, S=64, s=32, arc=2*pi){
mr <- transfoMatrix(p1,p2,p3)
tmesh0 <- createTorusMesh(mr$radius, r, S, s, arc)
transform3d(tmesh0, mr$matrix)
}
#### transformation matrix and radius ####
# plane passing by points p1, p2, p3 #
plane3pts <- function(p1,p2,p3){
xcoef = (p1[2]-p2[2])*(p2[3]-p3[3])-(p1[3]-p2[3])*(p2[2]-p3[2]);
ycoef = (p1[3]-p2[3])*(p2[1]-p3[1])-(p1[1]-p2[1])*(p2[3]-p3[3]);
zcoef = (p1[1]-p2[1])*(p2[2]-p3[2])-(p1[2]-p2[2])*(p2[1]-p3[1]);
offset = p1[1]*xcoef + p1[2]*ycoef + p1[3]*zcoef;
c(xcoef, ycoef, zcoef, offset)
}
# center, radius and normal of the circle passing by three points #
circleCenterAndRadius <- function(p1,p2,p3){
p12 <- (p1+p2)/2
p23 <- (p2+p3)/2
v12 <- p2-p1
v23 <- p3-p2
plane1 <- plane3pts(p1,p2,p3)
A <- rbind(plane1[1:3],v12,v23)
b <- c(plane1[4], sum(p12*v12), sum(p23*v23))
center <- c(solve(A) %*% b)
r <- sqrt(c(crossprod(p1-center)))
list(center=center, radius=r, normal=plane1[1:3])
}
# transformation matrix #
transfoMatrix <- function(p1,p2,p3){
crn <- circleCenterAndRadius(p1,p2,p3)
center <- crn$center
radius <- crn$radius
normal <- crn$normal
if(normal[1]==0 && normal[2]==0){
m <- rbind(cbind(diag(3), center), c(0,0,0,1))
}else{
measure <- sqrt(normal[1]^2+normal[2]^2+normal[3]^2)
normal <- normal/measure
s <- sqrt(normal[1]^2+normal[2]^2)
u <- c(normal[2]/s, -normal[1]/s, 0)
v <- geometry::extprod3d(normal, u)
m <- rbind(cbind(u, v, normal, center), c(0,0,0,1))
}
list(matrix=t(m), radius=radius)
}
#### torus as tmesh3d ####
# R: major radius
# r: minor radius
# S: number of segments for the major ring
# s: number of segments for the minor ring
# arc: the arc to draw
torusMesh <- function(R=1, r=0.25, S=64, s=32, arc=2*pi, addnormals=TRUE, ...){
vertices <- matrix(NA_real_, nrow = 3L, ncol = S*s)
indices1 <- indices2 <- matrix(NA_integer_, nrow = 3L, ncol = S*s)
u_ <- seq(0, arc, length.out = S+1)[-1]
v_ <- seq(0, 2*pi, length.out = s+1)[-1]
for(i in 1:S){
for(j in 1:s){
vertices[, (i-1)*s+j] <- c(
(R + r*cos(v_[j])) * cos(u_[i]),
(R + r*cos(v_[j])) * sin(u_[i]),
r * sin(v_[j])
)
}
}
if(addnormals){
Normals <- matrix(NA_real_, nrow = 3L, ncol = S*s)
for(i in 1L:S){
im1 <- ifelse(i==1L, S, i-1L)
ip1 <- ifelse(i==S, 1L, i+1L)
for(j in 1L:s){
jm1 <- ifelse(j==1L, s, j-1L)
jp1 <- ifelse(j==s, 1L, j+1L)
n1 <- rgl:::xprod(vertices[,(i-1)*s+jp1]-vertices[,(i-1)*s+j],
vertices[,(ip1-1)*s+j]-vertices[,(i-1)*s+j])
n2 <- -rgl:::xprod(vertices[,(i-1)*s+jm1]-vertices[,(i-1)*s+j],
vertices[,(ip1-1)*s+jm1]-vertices[,(i-1)*s+j])
n3 <- rgl:::xprod(vertices[,(im1-1)*s+j]-vertices[,(i-1)*s+j],
vertices[,(im1-1)*s+jp1]-vertices[,(i-1)*s+j])
n4 <- rgl:::xprod(vertices[,(ip1-1)*s+j]-vertices[,(i-1)*s+j],
vertices[,(ip1-1)*s+jm1]-vertices[,(i-1)*s+j])
n5 <- -rgl:::xprod(vertices[,(im1-1)*s+j]-vertices[,(i-1)*s+j],
vertices[,(i-1)*s+jm1]-vertices[,(i-1)*s+j])
n6 <- rgl:::xprod(vertices[,(im1-1)*s+jp1]-vertices[,(i-1)*s+j],
vertices[,(i-1)*s+jp1]-vertices[,(i-1)*s+j])
Normals[,(i-1)*s+j] <- c(n1+n2+n3+n4+n5+n6)/6
}
}
}
# triangles
s <- as.integer(s)
for(i in 1L:S){
ip1 <- ifelse(i==S, 1L, i+1L)
for(j in 1L:s){
jp1 <- ifelse(j==s, 1L, j+1L)
indices1[,(i-1)*s+j] <- c((i-1L)*s+j, (i-1L)*s+jp1, (ip1-1L)*s+j)
indices2[,(i-1)*s+j] <- c((i-1L)*s+jp1, (ip1-1L)*s+jp1, (ip1-1L)*s+j)
}
}
out <- list(
vb = rbind(vertices, 1),
it = cbind(indices1, indices2),
primitivetype = "triangle",
material = list(...),
normals = if(addnormals) rbind(Normals,1) else NULL
)
class(out) <- c("mesh3d", "shape3d")
out
}
#
torusMesh3points <- function(p1, p2, p3, r,
S=64, s=32, arc=2*pi, normals=TRUE, ...){
mr <- transfoMatrix(p1,p2,p3)
tmesh0 <- torusMesh(mr$radius, r, S, s, arc, normals, ...)
transform3d(tmesh0, mr$matrix)
}
# test
library(rgl)
p1 = c(1,-1,1)
p2 = c(2,1,2)
p3 = c(2,-2,-3)
tmesh <- torusMesh3points(p1, p2, p3, 0.5, color = "blue")
shade3d(tmesh)
triangles3d(rbind(p1,p2,p3), color="red")
#### transformation matrix and radius ####
# plane passing by points p1, p2, p3 #
plane3pts <- function(p1,p2,p3){
xcoef = (p1[2]-p2[2])*(p2[3]-p3[3])-(p1[3]-p2[3])*(p2[2]-p3[2]);
ycoef = (p1[3]-p2[3])*(p2[1]-p3[1])-(p1[1]-p2[1])*(p2[3]-p3[3]);
zcoef = (p1[1]-p2[1])*(p2[2]-p3[2])-(p1[2]-p2[2])*(p2[1]-p3[1]);
offset = p1[1]*xcoef + p1[2]*ycoef + p1[3]*zcoef;
c(xcoef, ycoef, zcoef, offset)
}
# center, radius and normal of the circle passing by three points #
circleCenterAndRadius <- function(p1,p2,p3){
p12 <- (p1+p2)/2
p23 <- (p2+p3)/2
v12 <- p2-p1
v23 <- p3-p2
plane1 <- plane3pts(p1,p2,p3)
A <- rbind(plane1[1:3],v12,v23)
b <- c(plane1[4], sum(p12*v12), sum(p23*v23))
center <- c(solve(A) %*% b)
r <- sqrt(c(crossprod(p1-center)))
list(center=center, radius=r, normal=plane1[1:3])
}
# transformation matrix #
transfoMatrix <- function(p1,p2,p3){
crn <- circleCenterAndRadius(p1,p2,p3)
center <- crn$center
radius <- crn$radius
normal <- crn$normal
if(normal[1]==0 && normal[2]==0){
m <- rbind(cbind(diag(3), center), c(0,0,0,1))
}else{
measure <- sqrt(normal[1]^2+normal[2]^2+normal[3]^2)
normal <- normal/measure
s <- sqrt(normal[1]^2+normal[2]^2)
u <- c(normal[2]/s, -normal[1]/s, 0)
v <- rgl:::xprod(normal, u)
m <- rbind(cbind(u, v, normal, center), c(0,0,0,1))
}
list(matrix=t(m), radius=radius)
}
#### torus as tmesh3d ####
# R: major radius
# r: minor radius
# S: number of subdivisions for the major ring
# s: number of subdivisions for the minor ring
torusMesh <- function(R=1, r=0.25, S=64, s=32, ...){
vertices <- matrix(NA_real_, nrow = 3L, ncol = S*s)
Normals <- matrix(NA_real_, nrow = 3L, ncol = S*s)
indices1 <- indices2 <- matrix(NA_integer_, nrow = 3L, ncol = S*s)
u_ <- seq(0, 2*pi, length.out = S+1)[-1]
v_ <- seq(0, 2*pi, length.out = s+1)[-1]
for(i in 1:S){
cos_ui <- cos(u_[i])
sin_ui <- sin(u_[i])
cx <- R * cos_ui
cy <- R * sin_ui
for(j in 1:s){
rcos_vj <- r*cos(v_[j])
n <- c(rcos_vj*cos_ui, rcos_vj*sin_ui, r*sin(v_[j]))
Normals[, (i-1)*s+j] <- -n
vertices[, (i-1)*s+j] <- c(cx,cy,0) + n
}
}
# triangles
s <- as.integer(s)
for(i in 1L:S){
ip1 <- ifelse(i==S, 1L, i+1L)
for(j in 1L:s){
jp1 <- ifelse(j==s, 1L, j+1L)
indices1[,(i-1)*s+j] <- c((i-1L)*s+j, (i-1L)*s+jp1, (ip1-1L)*s+j)
indices2[,(i-1)*s+j] <- c((i-1L)*s+jp1, (ip1-1L)*s+jp1, (ip1-1L)*s+j)
}
}
out <- list(
vb = rbind(vertices, 1),
it = cbind(indices1, indices2),
primitivetype = "triangle",
material = list(...),
normals = rbind(Normals,1)
)
class(out) <- c("mesh3d", "shape3d")
out
}
#
torusMesh3points <- function(p1, p2, p3, r,
S=64, s=32, normals=TRUE, ...){
mr <- transfoMatrix(p1,p2,p3)
tmesh0 <- torusMesh(mr$radius, r, S, s, normals, ...)
transform3d(tmesh0, mr$matrix)
}
# # test
# library(rgl)
# p1 = c(1,-1,1)
# p2 = c(2,1,2)
# p3 = c(2,-2,-3)
# tmesh <- torusMesh3points(p1, p2, p3, 0.5, color = "blue")
# shade3d(tmesh)
# triangles3d(rbind(p1,p2,p3), color="red")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment