Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- setwd("~/Work/R/rgl/Meshes")
- library(rgl)
- # knot curve
- knot <- function(t, p, q){
- r <- cos(q*t)+2
- c(
- r * cos(p*t),
- r * sin(p*t),
- -sin(q*t)
- )
- }
- # derivative (tangent)
- dknot <- function(t, p, q){
- v <- c(
- -q*sin(q*t)*cos(p*t) - p*sin(p*t)*(cos(q*t)+2),
- -q*sin(q*t)*sin(p*t) + p*cos(p*t)*(cos(q*t)+2),
- -q*cos(q*t)
- )
- v / sqrt(c(crossprod(v)))
- }
- # second derivative (normal)
- ddknot <- function(t, p, q){
- v <- c(
- -q*(q*cos(q*t)*cos(p*t)-p*sin(p*t)*sin(q*t)) -
- p*(p*cos(p*t)*cos(q*t+2)-q*sin(q*t)*sin(p*t)),
- -q*(q*cos(q*t)*sin(p*t)+p*cos(p*t)*sin(q*t)) +
- p*(-p*sin(p*t)*cos(q*t+2)-q*sin(q*t)*cos(p*t)),
- q*q*sin(q*t)
- )
- v / sqrt(c(crossprod(v)))
- }
- # binormal
- bnrml <- function(t, p, q){
- v <- rgl:::xprod(dknot(t, p, q), ddknot(t, p, q))
- v / sqrt(c(crossprod(v)))
- }
- # mesh: twisted tubular astroknot
- scos <- function(x,alpha) sign(cos(x)) * abs(cos(x))^alpha
- ssin <- function(x,alpha) sign(sin(x)) * abs(sin(x))^alpha
- TubularKnotMesh <- function(p, q, a, nu, nv, alpha=1, twists=2){
- vs <- matrix(NA_real_, nrow=3, ncol=nu*nv)
- colors <- matrix(NA_character_, nrow = 3L, ncol = nu*nv)
- u_ <- seq(0, 2*pi, length.out = nu+1)[-1]
- v_ <- seq(0, 2*pi, length.out = nv+1)[-1]
- for(i in 1:nu){
- u <- u_[i]
- for(j in 1:nv){
- v <- v_[j]
- h <- knot(u, p, q)
- vs[,(i-1)*nv+j] <-
- h +
- a*(scos(v,alpha) *
- (cos(twists*u)*ddknot(u, p, q) +
- sin(twists*u)*bnrml(u, p, q)) +
- ssin(v,alpha) *
- (-sin(twists*u)*ddknot(u, p, q) +
- cos(twists*u)*bnrml(u, p, q)))
- colors[,(i-1)*nv+j] <-
- ifelse(findInterval(v, pi*(1/16+seq(0,2,len=17))) %% 2 == 0,
- "#440154FF", "#FDE725FF")
- }
- }
- tris1 <- matrix(NA_integer_, nrow=3, ncol=nu*nv)
- tris2 <- matrix(NA_integer_, nrow=3, ncol=nu*nv)
- nv <- as.integer(nv)
- for(i in 1L:nu){
- ip1 <- ifelse(i==nu, 1L, i+1L)
- for(j in 1L:nv){
- jp1 <- ifelse(j==nv, 1L, j+1L)
- tris1[,(i-1)*nv+j] <- c((i-1L)*nv+j,(i-1L)*nv+jp1, (ip1-1L)*nv+j)
- tris2[,(i-1)*nv+j] <- c((i-1L)*nv+jp1,(ip1-1L)*nv+jp1,(ip1-1L)*nv+j)
- }
- }
- out <- tmesh3d(
- vertices = vs,
- indices = cbind(tris1, tris2),
- homogeneous = FALSE,
- material = list(color = colors)
- )
- addNormals(out)
- }
- m <- TubularKnotMesh(p=3, q=5, a=0.5, 10*60, 60, alpha=1, twists=3)
- open3d(windowRect = c(50,50,550,550))
- bg3d(rgb(54,57,64, maxColorValue = 255))
- view3d(0, 0, zoom=0.75)
- shade3d(m)
- # animation ####
- movie3d(spin3d(axis = c(0, 0, 1), rpm = 60),
- duration = 1, fps = 60,
- movie = "anim", dir = ".",
- convert = "convert -dispose previous -loop 0 -delay 1x%d %s*.png %s.%s",
- startTime = 1/60)
Add Comment
Please, Sign In to add comment