I previously asked a question on creating 3D dendrograms in rgl window in R and @user2554330 provided the solution here. I then modified @user2554330's code to add 3D meshes to the 3D dendrogram.
Fully reproducible code:
required<-c("rgl", "ggdendro", "Morpho")
if(any(!required%in%installed.packages()[,1])){
install.packages(required[which(!required%in%installed.packages()[,1])])
}
sapply(required,require,character.only=TRUE)
# Create mesh
mshp <- matrix(c(
-1.0, 0.0, 0.0,
1.0, 0.0, 0.0,
0.0, -1.0, 0.0,
0.0, 1.0, 0.0,
0.0, 0.0, -1.0,
0.0, 0.0, 1.0
), ncol = 3, byrow = T)
itFull <- matrix(c(
1,5,3,
1,3,6,
1,4,5,
1,6,4,
2,3,5,
2,6,3,
2,5,4,
2,4,6
), ncol = 3, byrow = T)
# Center the mesh.
mshp <- scale(mshp, center = T, scale = F)
lm2mesh <- function(vb, it, addNormals = T) {
library(rgl)
vb <- t(vb)
vb <- rbind(vb, 1)
rownames(vb) <- c("xpts", "ypts", "zpts", "")
it_mat <- t(as.matrix(it))
rownames(it_mat) <- NULL
vertices <- c(vb)
indices <- c(it_mat)
meshExp <- tmesh3d(vertices = vertices, indices = indices, homogeneous = TRUE,
material = NULL, normals = NULL, texcoords = NULL)
if (addNormals == T) {
meshExp <- addNormals(meshExp)
}
return(meshExp)
}
mMesh <- lm2mesh(mshp, itFull)
# Scale mesh to unit size.
mMesh <- scalemesh(mMesh, 1/cSize(mMesh), center = "none")
################################################################################
# Construct dendrogram in 3D
################################################################################
hclust <- list()
hclust$merge <- matrix(c(-1, -2,
-3, 1,
-4, 2,
-5, 3,
-6, -7,
-8, -9,
5, 6,
-12, -13,
-11, 8,
-10, 9,
7, 10,
4, 11), ncol = 2, byrow = T)
hclust$height <- c(1-0.8702480, 1-0.8359299, 1-0.7754834, 1-0.5989334,
1-0.8289807, 1-0.8142845, 1-0.6319150, 1-0.7453669,
1-0.5331113, 1-0.5331113, 1-0.2472344, 1-0.2472344)
hclust$order <- 1:13 # order of leaves(trivial if hand-entered)
hclust$labels <- LETTERS[1:13] # labels of leaves
class(hclust) <- "hclust" # make it an hclust object
plot(hclust)
# Convert to a dendrogram object.
hclustd <- as.dendrogram(hclust)
dend_data <- dendro_data(hclustd, type = "rectangle")
nodes <- dend_data$segments
# Set the gap between the ends of the tree
gap <- 0
# Set the offset from the center.
offset <- 0
radius <- with(nodes, max(c(y, yend)) + offset)
circ <- with(nodes, max(c(x, xend)) + gap)
# Convert to polar coordinates
nodes$theta <- with(nodes, 2*pi*x/circ)
nodes$thetaend <- with(nodes, 2*pi*xend/circ)
nodes$r <- with(nodes, (radius - y)/radius)
nodes$rend <- with(nodes, (radius - yend)/radius)
# Extract the horizontal and vertical segments
horiz <- subset(nodes, y == yend)
vert <- subset(nodes, x == xend)
library(rgl)
open3d(userMatrix = diag(4))
# Draw the vertical segments, which are still segments
x <- with(vert, as.numeric(rbind(r*cos(theta), rend*cos(theta))))
y <- with(vert, as.numeric(rbind(r*sin(theta), rend*sin(theta))))
segments3d(x, y, z = -0.1)
# Draw the horizontal segments, which are now arcs. Zero
# radius arcs are dropped
horiz <- subset(horiz, r > 0)
# For row 17, x=xend and y=yend, so it is a point. This leads to problems when
# using arc3d.
with(horiz[c(1:16, 18:21),], arc3d(from = cbind(r*cos(theta), r*sin(theta), -0.1),
to = cbind(r*cos(thetaend), r*sin(thetaend), -0.1),
center = c(0, 0, -0.1)))
# Draw the labels
labels <- dend_data$labels
labels$theta <- with(labels, 2*pi*x/circ)
# Add a bit to the y so the label doesn't overlap the segment
labels$r <- with(labels, (radius - y)/radius + 0.1)
with(labels, text3d(r*cos(theta), r*sin(theta), 0, label))
################################################################################
# Plot facial mesh
################################################################################
# Function that adds mesh to desired location.
addMesh <- function(useMesh, x_center, y_center, z_center, scl = c(0.5, 0.5, 0.5),
meshColor = c("vertices", "faces"), specular = "#202020",
alpha = 0.9,
open3d = F) {
# Scale the useMesh object to a desired size (if needed), and translate it to
# the desired position.
useMesh$vb[1:3,] <- useMesh$vb[1:3,] * scl[row(useMesh$vb[1:3,])]
useMesh$vb[1,] = useMesh$vb[1,] + x_center
useMesh$vb[2,] = useMesh$vb[2,] + y_center
useMesh$vb[3,] = useMesh$vb[3,] + z_center
newMesh <- shade3d(useMesh, meshColor = meshColor, alpha = alpha,
color = "lightgrey", specular = specular, open3d = open3d)
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Add meshes to dendrogram one by one.
# Module #6 (The first module when Sall[,1] is sorted in ascending order).
newMesh6 <- addMesh(mMesh,
x_center = (vert$rend[2] + 0.2) * cos(vert$thetaend[2]),
y_center = (vert$rend[2] + 0.2) * sin(vert$thetaend[2]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #9
newMesh9 <- addMesh(mMesh,
x_center = (vert$rend[19] + 0.2) * cos(vert$thetaend[19]),
y_center = (vert$rend[19] + 0.2) * sin(vert$thetaend[19]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #10
newMesh10 <- addMesh(mMesh,
x_center = (vert$rend[21] + 0.2) * cos(vert$thetaend[21]),
y_center = (vert$rend[21] + 0.2) * sin(vert$thetaend[21]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #12
newMesh12 <- addMesh(mMesh,
x_center = (vert$rend[4] + 0.2) * cos(vert$thetaend[4]),
y_center = (vert$rend[4] + 0.2) * sin(vert$thetaend[4]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #14
newMesh14 <- addMesh(mMesh,
x_center = (vert$rend[13] + 0.2) * cos(vert$thetaend[13]),
y_center = (vert$rend[13] + 0.2) * sin(vert$thetaend[13]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #15
newMesh15 <- addMesh(mMesh,
x_center = (vert$rend[14] + 0.2) * cos(vert$thetaend[14]),
y_center = (vert$rend[14] + 0.2) * sin(vert$thetaend[14]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #16
newMesh16 <- addMesh(mMesh,
x_center = (vert$rend[16] + 0.2) * cos(vert$thetaend[16]),
y_center = (vert$rend[16] + 0.2) * sin(vert$thetaend[16]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #17
newMesh17 <- addMesh(mMesh,
x_center = (vert$rend[17] + 0.2) * cos(vert$thetaend[17]),
y_center = (vert$rend[17] + 0.2) * sin(vert$thetaend[17]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #18
newMesh18 <- addMesh(mMesh,
x_center = (vert$rend[23] + 0.2) * cos(vert$thetaend[23]),
y_center = (vert$rend[23] + 0.2) * sin(vert$thetaend[23]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #19
newMesh19 <- addMesh(mMesh,
x_center = (vert$rend[24] + 0.2) * cos(vert$thetaend[24]),
y_center = (vert$rend[24] + 0.2) * sin(vert$thetaend[24]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #20
newMesh20 <- addMesh(mMesh,
x_center = (vert$rend[6] + 0.2) * cos(vert$thetaend[6]),
y_center = (vert$rend[6] + 0.2) * sin(vert$thetaend[6]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #22
newMesh22 <- addMesh(mMesh,
x_center = (vert$rend[8] + 0.2) * cos(vert$thetaend[8]),
y_center = (vert$rend[8] + 0.2) * sin(vert$thetaend[8]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #23
newMesh23 <- addMesh(mMesh,
x_center = (vert$rend[9] + 0.2) * cos(vert$thetaend[9]),
y_center = (vert$rend[9] + 0.2) * sin(vert$thetaend[9]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Outermost circle completed.
# Module #21
newMesh21 <- addMesh(mMesh,
x_center = (vert$rend[7]) * cos(vert$thetaend[7]),
y_center = (vert$rend[7]) * sin(vert$thetaend[7]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #13
newMesh13 <- addMesh(mMesh,
x_center = (vert$rend[5]) * cos(vert$thetaend[5] - 0.1),
y_center = (vert$rend[5]) * sin(vert$thetaend[5] - 0.1),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #5
newMesh5 <- addMesh(mMesh,
x_center = (vert$rend[3]) * cos(vert$thetaend[3] - 0.15),
y_center = (vert$rend[3]) * sin(vert$thetaend[3] - 0.15),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #2
newMesh2 <- addMesh(mMesh,
x_center = (vert$rend[1] - 0.1) * cos(vert$thetaend[1]),
y_center = (vert$rend[1] - 0.1) * sin(vert$thetaend[1]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #7
newMesh7 <- addMesh(mMesh,
x_center = (vert$rend[12]) * cos(vert$thetaend[12]),
y_center = (vert$rend[12]) * sin(vert$thetaend[12]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #8
newMesh8 <- addMesh(mMesh,
x_center = (vert$rend[15]) * cos(vert$thetaend[15]),
y_center = (vert$rend[15]) * sin(vert$thetaend[15]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #3
newMesh3 <- addMesh(mMesh,
x_center = (vert$rend[11]) * cos(vert$thetaend[11]),
y_center = (vert$rend[11]) * sin(vert$thetaend[11]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #11
newMesh11 <- addMesh(mMesh,
x_center = (vert$rend[22]) * cos(vert$thetaend[22]),
y_center = (vert$rend[22]) * sin(vert$thetaend[22]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Module #4
newMesh4 <- addMesh(mMesh,
x_center = (vert$rend[18]) * cos(vert$thetaend[18]),
y_center = (vert$rend[18]) * sin(vert$thetaend[18]),
z_center = 0,
meshColor = "faces",
open3d = F)
# Done!
The angles of the tetrahedron is slightly different. I want all meshes to be in the same orientation.
###################################################
@user2554330 pointed out in comment that this is due to parallex issue. I added par3d(FOV=0)
to the last line of code and now all meshes are oriented identically. My issue is solved.
Although par3d(FOV=0)
perfectly solved my issue, each mesh cannot be rotated individually. I wish to further modify the figure to allow for rotation of each mesh individually.
A possible solution would be to create subscenes at locations surrounding the tetrahedrons and specify the same orientation for all meshes in their respective subscenes. @user2554330 provided the method to create subscene in rgl window in another of my question. However, I do not know how to specify newviewport
and vert$theta
and vert$thetaend
for each mesh subscene.