Skip to content

Instantly share code, notes, and snippets.

@stla
Last active Nov 27, 2018
Embed
What would you like to do?
Spherical triangle with rgl - https://i.imgur.com/05qwdir.gifv
library(rgl)
library(abind) # to use the abind function
library(plyr) # to use the alply function
# normalize a vector
normalize <- function(vec){
vecnorm <- sqrt(sum(vec^2))
vec/vecnorm
}
# split a triangle into 4 triangles
# vertices of the triangle are assumed to lie on the unit sphere
splitTriangle1 <- function(triangle){
a <- triangle[1,]
b <- triangle[2,]
c <- triangle[3,]
mab <- normalize((a+b)/2)
mac <- normalize((a+c)/2)
mbc <- normalize((b+c)/2)
abind(
rbind(a,mab,mac),
rbind(b,mbc,mab),
rbind(c,mac,mbc),
rbind(mab,mbc,mac),
along = 3
)
}
# example of splitTriangle1
triangle <- rbind(
c(0,0,1),
c(0,1,0),
c(1,0,0)
)
splitTriangle1(triangle)
# iterate the splitTriangle1 function
splitTriangle <- function(triangle, n){
triangles <- splitTriangle1(triangle)
for(i in seq_len(n-1)){
triangles <- abind(alply(triangles,3, splitTriangle1), along=3)
}
unname(triangles)
}
# example of splitTriangle
triangle <- rbind(
c(1,0,0),
c(0,1,0),
c(0,0,1)
)
triangles <- splitTriangle(triangle, 3)
for(i in seq_len(dim(triangles)[3])){
triangles3d(triangles[,,i], color="green")
}
# main function: draw a spherical triangle
drawSphericalTriangle <- function(triangle, n=5, col="green"){
triangles <- splitTriangle(triangle, n)
for(i in seq_len(dim(triangles)[3])){
triangles3d(triangles[,,i], color=col)
}
}
#### example: triangulated sphere ####
randomOnSphere <- function(n){
theta <- runif(n, 0, 2*pi)
phi <- runif(n, 0, pi)
cbind(cos(theta)*sin(phi), sin(theta)*sin(phi), cos(phi))
}
points <- randomOnSphere(50)
library(geometry)
hull <- convhulln(points, options="Qt")
triangles <- alply(hull, 1,
function(ijk){
rbind(points[ijk[1],],points[ijk[2],],points[ijk[3],])
})
colors <- randomcoloR::distinctColorPalette(length(triangles))
lapply(seq_along(triangles),
function(i) drawSphericalTriangle(triangles[[i]], n=5, col=colors[i]))
library(rgl)
library(abind) # to use the abind function
library(plyr) # to use the alply function
# normalize a vector
normalize <- function(vec){
vecnorm <- sqrt(sum(vec^2))
vec/vecnorm
}
# split a triangle into 4 triangles
# vertices of the triangle are assumed to lie on the unit sphere
splitTriangle1 <- function(triangle){
a <- triangle[1,]
b <- triangle[2,]
c <- triangle[3,]
mab <- normalize((a+b)/2)
mac <- normalize((a+c)/2)
mbc <- normalize((b+c)/2)
abind(
rbind(a,mab,mac),
rbind(b,mbc,mab),
rbind(c,mac,mbc),
rbind(mab,mbc,mac),
along = 3
)
}
# example of splitTriangle1
triangle <- rbind(
c(0,0,1),
c(0,1,0),
c(1,0,0)
)
splitTriangle1(triangle)
# iterate the splitTriangle1 function to do the mesh
sphtriMesh <- function(triangle, n){
normal <- rgl:::xprod(triangle[2,]-triangle[1,],triangle[3,]-triangle[1,])
center <- colMeans(triangle)
dot <- crossprod(normal, center)
if(dot < 0){
triangle <- triangle[c(1,3,2),]
}
triangles <- splitTriangle1(triangle)
for(i in seq_len(n-1)){
triangles <- abind(alply(triangles,3, splitTriangle1), along=3)
}
vertices <- do.call(cbind, alply(triangles, 3, t))
tris <- matrix(1:(3*4^n), nrow=3)
out <- list(
vb = rbind(vertices,1),
it = tris,
primitivetype = "triangle",
material = list(),
normals = rbind(vertices,1)
)
class(out) <- c("mesh3d", "shape3d")
out
}
# example
triangle <- rbind(
c(1,0,0),
c(0,1,0),
c(0,0,1)
)
mesh <- sphtriMesh(triangle, 4)
shade3d(mesh, color="red")
triangle <- rbind(
c(1,0,0),
c(0,0,1),
c(0,1,0)
)
mesh <- sphtriMesh(triangle, 4)
shade3d(mesh, color="red")
#### example: triangulated sphere ####
points <- uniformly::runif_on_sphere(50, 3)
library(geometry)
hull <- convhulln(points, options="Qt")
triangles <- alply(hull, 1,
function(ijk){
rbind(points[ijk[1],],points[ijk[2],],points[ijk[3],])
})
colors <- randomcoloR::distinctColorPalette(length(triangles))
rgl.material(back="culled")
invisible(lapply(seq_along(triangles),
function(i){
mesh <- sphtriMesh(triangles[[i]], n=3)
shade3d(mesh, color=colors[i])
}))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment