• API
• FAQ
• Tools
• Archive
daily pastebin goal
41%
SHARE
TWEET

# Untitled

a guest Dec 12th, 2018 58 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.   )
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)