daily pastebin goal
68%
SHARE
TWEET

Untitled

a guest Dec 12th, 2018 59 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. setwd("~/Work/R/rgl/Meshes")
  2. library(rgl)
  3.  
  4. # knot curve
  5. knot <- function(t, p, q){
  6.   r <- cos(q*t)+2
  7.   c(
  8.     r * cos(p*t),
  9.     r * sin(p*t),
  10.     -sin(q*t)
  11.   )
  12. }
  13. # derivative (tangent)
  14. dknot <- function(t, p, q){
  15.   v <- c(
  16.     -q*sin(q*t)*cos(p*t) - p*sin(p*t)*(cos(q*t)+2),
  17.     -q*sin(q*t)*sin(p*t) + p*cos(p*t)*(cos(q*t)+2),
  18.     -q*cos(q*t)
  19.   )
  20.   v / sqrt(c(crossprod(v)))
  21. }
  22. # second derivative (normal)
  23. ddknot <- function(t, p, q){
  24.   v <- c(
  25.     -q*(q*cos(q*t)*cos(p*t)-p*sin(p*t)*sin(q*t)) -
  26.       p*(p*cos(p*t)*cos(q*t+2)-q*sin(q*t)*sin(p*t)),
  27.     -q*(q*cos(q*t)*sin(p*t)+p*cos(p*t)*sin(q*t)) +
  28.       p*(-p*sin(p*t)*cos(q*t+2)-q*sin(q*t)*cos(p*t)),
  29.     q*q*sin(q*t)
  30.   )
  31.   v / sqrt(c(crossprod(v)))
  32. }
  33. # binormal
  34. bnrml <- function(t, p, q){
  35.   v <- rgl:::xprod(dknot(t, p, q), ddknot(t, p, q))
  36.   v / sqrt(c(crossprod(v)))
  37. }
  38.  
  39. # mesh: twisted tubular astroknot
  40. scos <- function(x,alpha) sign(cos(x)) * abs(cos(x))^alpha
  41. ssin <- function(x,alpha) sign(sin(x)) * abs(sin(x))^alpha
  42. TubularKnotMesh <- function(p, q, a, nu, nv, alpha=1, twists=2){
  43.   vs <- matrix(NA_real_, nrow=3, ncol=nu*nv)
  44.   colors <- matrix(NA_character_, nrow = 3L, ncol = nu*nv)
  45.   u_ <- seq(0, 2*pi, length.out = nu+1)[-1]
  46.   v_ <- seq(0, 2*pi, length.out = nv+1)[-1]
  47.   for(i in 1:nu){
  48.     u <- u_[i]
  49.     for(j in 1:nv){
  50.       v <- v_[j]
  51.       h <- knot(u, p, q)
  52.       vs[,(i-1)*nv+j] <-
  53.         h +
  54.         a*(scos(v,alpha) *
  55.              (cos(twists*u)*ddknot(u, p, q) +
  56.                 sin(twists*u)*bnrml(u, p, q)) +
  57.            ssin(v,alpha) *
  58.              (-sin(twists*u)*ddknot(u, p, q) +
  59.                 cos(twists*u)*bnrml(u, p, q)))
  60.       colors[,(i-1)*nv+j] <-
  61.         ifelse(findInterval(v, pi*(1/16+seq(0,2,len=17))) %% 2 == 0,
  62.                "#440154FF", "#FDE725FF")
  63.     }
  64.   }
  65.   tris1 <- matrix(NA_integer_, nrow=3, ncol=nu*nv)
  66.   tris2 <- matrix(NA_integer_, nrow=3, ncol=nu*nv)
  67.   nv <- as.integer(nv)
  68.   for(i in 1L:nu){
  69.     ip1 <- ifelse(i==nu, 1L, i+1L)
  70.     for(j in 1L:nv){
  71.       jp1 <- ifelse(j==nv, 1L, j+1L)
  72.       tris1[,(i-1)*nv+j] <- c((i-1L)*nv+j,(i-1L)*nv+jp1, (ip1-1L)*nv+j)  
  73.       tris2[,(i-1)*nv+j] <- c((i-1L)*nv+jp1,(ip1-1L)*nv+jp1,(ip1-1L)*nv+j)  
  74.     }
  75.   }
  76.   out <- tmesh3d(
  77.     vertices = vs,
  78.     indices = cbind(tris1, tris2),
  79.     homogeneous = FALSE,
  80.     material = list(color = colors)
  81.   )
  82.   addNormals(out)
  83. }
  84.  
  85. m <- TubularKnotMesh(p=3, q=5, a=0.5, 10*60, 60, alpha=1, twists=3)
  86. open3d(windowRect = c(50,50,550,550))
  87. bg3d(rgb(54,57,64, maxColorValue = 255))
  88. view3d(0, 0, zoom=0.75)
  89. shade3d(m)
  90.  
  91. # animation ####
  92. movie3d(spin3d(axis = c(0, 0, 1), rpm = 60),
  93.         duration = 1, fps = 60,
  94.         movie = "anim", dir = ".",
  95.         convert = "convert -dispose previous -loop 0 -delay 1x%d %s*.png %s.%s",
  96.         startTime = 1/60)
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top