Torus passing by three points in rgl

TorusPassingByThreePoints.R

# 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)
}

TorusPassingByThreePoints2.R

#### 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")

TorusPassingByThreePoints3.R

#### 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")