Drawing a torus with rgl

Posted on May 29, 2018 by Stéphane Laurent
Tags: R, graphics, rgl, geometry

The following function creates a torus as a mesh3d object, which can be drawn with the rgl package.

# R: major 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) {
for (i in 0:S) {
u <- i / S * arc
v <- j / s * 2*pi
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
c <- (S + 1) * (j - 1) + i + 1
d <- (S + 1) * j + i + 1
indices <- cbind(indices, c(a,b,d), c(b,c,d))
}
}
tmesh3d(rbind(vertices,1), indices)
}
torus <- createTorusMesh(R=2, r=1)
shade3d(torus, color="blue")

Denote by $$O$$ the center of the torus and by $$R$$ its major radius. The equator of this torus (circle of center $$O$$ and radius $$R$$) lies in the $$xy$$-plane (the plane $$z=0$$).

Now we are going to give a method to draw a torus such that its equator passes by three given points. This will be done with the help of the transform3d function of the rgl package. Thus, the key element to construct is the transformation matrix.

Firstly, we write a function which returns the center and the radius of a circle passing by three given points in the 3D space. This is the function circleCenterAndRadius below. The plane3pts function is an helper function which returns the equation of the plane passing by three points.

# 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 #
p12 <- (p1+p2)/2
p23 <- (p2+p3)/2
v12 <- p2-p1
v23 <- p3-p2
plane <- plane3pts(p1,p2,p3)
A <- rbind(plane[1:3],v12,v23)
b <- c(plane[4], sum(p12*v12), sum(p23*v23))
center <- c(solve(A) %*% b)
r <- sqrt(c(crossprod(p1-center)))
}
Note that circleCenterAndRadius actually returns not only the center and the radius, but also a vector perpendicular to the plane of the circle.
transfoMatrix <- function(p1,p2,p3){
center <- crn$center radius <- crn$radius
normal <- crn$normal measure <- sqrt(normal[1]^2+normal[2]^2+normal[3]^2) normal <- normal/measure s <- sqrt(normal[1]^2+normal[2]^2) # TODO: case s=0 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) } Our function transfoMatrix also returns a radius, the major radius of the desired torus. Let us test it now. library(rgl) # our three points p1 = c(1,-1,1) p2 = c(2,1,2) p3 = c(2,-2,-3) # get the transformation matrix and the radius mr <- transfoMatrix(p1,p2,p3) tmesh0 <- createTorusMesh(R=mr$radius, r=0.1)
tmesh <- transform3d(tmesh0, mr$matrix) # draw the torus shade3d(tmesh, color="blue") # and draw a triangle to check triangles3d(rbind(p1,p2,p3), color="red") It works. So let’s play with it now. tetrahedron <- tetrahedron3d()[c("vb","it")] vertices <- tetrahedron$vb[1:3,]
triangles <- tetrahedron$it # p1 <- vertices[,triangles[1,1]] p2 <- vertices[,triangles[2,1]] p3 <- vertices[,triangles[3,1]] mr <- transfoMatrix(p1,p2,p3) tmesh <- transform3d(createTorusMesh(R=mr$radius, r=0.1), mr$matrix) shade3d(tmesh, color="blue") # p1 <- vertices[,triangles[1,2]] p2 <- vertices[,triangles[2,2]] p3 <- vertices[,triangles[3,2]] mr <- transfoMatrix(p1,p2,p3) tmesh <- transform3d(createTorusMesh(R=mr$radius, r=0.1), mr$matrix) shade3d(tmesh, color="blue") # p1 <- vertices[,triangles[1,3]] p2 <- vertices[,triangles[2,3]] p3 <- vertices[,triangles[3,3]] mr <- transfoMatrix(p1,p2,p3) tmesh <- transform3d(createTorusMesh(R=mr$radius, r=0.1), mr$matrix) shade3d(tmesh, color="blue") # p1 <- vertices[,triangles[1,4]] p2 <- vertices[,triangles[2,4]] p3 <- vertices[,triangles[3,4]] mr <- transfoMatrix(p1,p2,p3) tmesh <- transform3d(createTorusMesh(R=mr$radius, r=0.1), mr\$matrix)
shade3d(tmesh, color="blue")