Guest User

Untitled

a guest
Dec 12th, 2018
107
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.68 KB | None | 0 0
  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)
Add Comment
Please, Sign In to add comment