在 R

Question

我之前问过一个关于在 R 的 rgl 窗口中创建 3D 树状图的问题，@user2554330 提供了解决方案 here。然后我修改了@user2554330 的代码，将 3D 面部网格添加到 3D 树状图中。

下面是代码。可以在这里访问数据。 （面部网格来自公共资源库，不涉及隐私问题）

# Load libraries
required<-c("rgl", "ggdendro", "Morpho", "rmatio", "colorspace", "readxl", "data.table")
if(any(!required%in%installed.packages()[,1])){
    install.packages(required[which(!required%in%installed.packages()[,1])])
}
sapply(required,require,character.only=TRUE)

# Put all 6 downloaded files into a folder named 3DTree (remember to modify your path to the 3DTree folder)
mypath <- "<YOUR PATH HERE>\\3DTree"

# Load custom functions, which are used to create facial mesh, superimpose 
# colors, and add to specified position in rgl window.
source(paste0(mypath, "\\CustomeFunction.R"))

################################################################################
# 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))


################################################################################
# Add facial mesh
################################################################################
# Load facial segmentation information
Sall <- read.mat(paste0(mypath, "\\Sall.mat"))
Sall <- Sall$Sall
# Column 3 of Tree contains values to be mapped to color
ScTree <- read.mat(paste0(mypath, "\\ScTree.mat"))
Tree <- ScTree$Tree

# Pick color
palPurple <- sequential_hcl(n = 40, h = 270, c = c(70, NA, NA), l = c(25, 95), power = 1.2, register = )
palPurple <- rev(colorRampPalette(palPurple)(100))
palRed <- sequential_hcl(n = 40, h = 10, c = c(85, NA, NA), l = c(25, 95), power = 1.3, register = )
palRed <- rev(colorRampPalette(palRed)(100))
palGreen <- sequential_hcl(n = 40, h = 135, c = c(45, NA, NA), l = c(35, 95), power = 1.3, register = )
palGreen <- rev(colorRampPalette(palGreen)(100))

# Load full facial mesh information.
mshp <- read_excel(paste0(mypath, "\\Template.xlsx"), sheet = "Sheet1", col_names = F)
itFull <- read_excel(paste0(mypath, "\\it.xlsx"), sheet = "Sheet1", col_names = F)

# Rotate loaded facial mesh 30 degrees to the left.
# rotMat_y rotates faces around the y-axis, i.e., within the x-z plane.
theta_y <- -30 * pi/180
rotMat_y <- matrix(c( cos(theta_y), 0,  sin(theta_y),
                      0,            1,  0,
                     -sin(theta_y), 0,  cos(theta_y)), 
                   byrow = T, ncol = 3)
mshp <- as.matrix(mshp) %*% rotMat_y

# Center the mesh.
mshp <- scale(mshp, center = T, scale = F)
mMesh <- lm2mesh(mshp, itFull)
# Scale mesh to unit size.
mMesh <- scalemesh(mMesh, 1/cSize(mMesh), center = "none")

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Add meshes to dendrogram one by one.
# Module #6 (The first module when Sall[,1] is sorted in ascending order).
colHighlight6 <- palRed[round(Tree[which(Tree[,2] == 6), 3]*100, 0)]
palUse6 <- rep("#e6e6e6", length(unique(Sall[,1])))
palUse6[which(sort(unique(Sall[,1])) == 6)] <- colHighlight6
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,
                   scl = c(12, 12, 12),
                   locDiff = locDiffMirror(match(Sall[,1], sort(unique(Sall[,1])))),
                   from = 1, to = length(unique(Sall[,1])), 
                   pal = palUse6, 
                   meshColor = "faces", 
                   open3d = F)

# Module #9
colHighlight9 <- palGreen[round(Tree[which(Tree[,2] == 9), 3]*100, 0)]
palUse9 <- rep("#e6e6e6", length(unique(Sall[,1])))
palUse9[which(sort(unique(Sall[,1])) == 9)] <- colHighlight9
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,
                   scl = c(12, 12, 12),
                   locDiff = match(Sall[,1], sort(unique(Sall[,1]))),
                   from = 1, to = length(unique(Sall[,1])), 
                   pal = palUse9, 
                   meshColor = "faces", 
                   open3d = F)

# Module #10
colHighlight10 <- palGreen[round(Tree[which(Tree[,2] == 10), 3]*100, 0)]
palUse10 <- rep("#e6e6e6", length(unique(Sall[,1])))
palUse10[which(sort(unique(Sall[,1])) == 10)] <- colHighlight10
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,
                   scl = c(12, 12, 12),
                   locDiff = match(Sall[,1], sort(unique(Sall[,1]))),
                   from = 1, to = length(unique(Sall[,1])), 
                   pal = palUse10, 
                   meshColor = "faces", 
                   open3d = F)

# Module #12
colHighlight12 <- palRed[round(Tree[which(Tree[,2] == 12), 3]*100, 0)]
palUse12 <- rep("#e6e6e6", length(unique(Sall[,1])))
palUse12[which(sort(unique(Sall[,1])) == 12)] <- colHighlight12
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,
                   scl = c(12, 12, 12),
                   locDiff = match(Sall[,1], sort(unique(Sall[,1]))),
                   from = 1, to = length(unique(Sall[,1])), 
                   pal = palUse12, 
                   meshColor = "faces", 
                   open3d = F)

# Module #14
colHighlight14 <- palPurple[round(Tree[which(Tree[,2] == 14), 3]*100, 0)]
palUse14 <- rep("#F0F0F0", length(unique(Sall[,1])))
palUse14[which(sort(unique(Sall[,1])) == 14)] <- colHighlight14
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,
                   scl = c(12, 12, 12),
                   locDiff = match(Sall[,1], sort(unique(Sall[,1]))),
                   from = 1, to = length(unique(Sall[,1])), 
                   pal = palUse14, 
                   meshColor = "faces", 
                   open3d = F)

# Module #15
colHighlight15 <- palPurple[round(Tree[which(Tree[,2] == 15), 3]*100, 0)]
palUse15 <- rep("#F0F0F0", length(unique(Sall[,1])))
palUse15[which(sort(unique(Sall[,1])) == 15)] <- colHighlight15
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,
                   scl = c(12, 12, 12),
                   locDiff = match(Sall[,1], sort(unique(Sall[,1]))),
                   from = 1, to = length(unique(Sall[,1])), 
                   pal = palUse15, 
                   meshColor = "faces", 
                   open3d = F)

# Module #16
colHighlight16 <- palPurple[round(Tree[which(Tree[,2] == 16), 3]*100, 0)]
palUse16 <- rep("#F0F0F0", length(unique(Sall[,1])))
palUse16[which(sort(unique(Sall[,1])) == 16)] <- colHighlight16
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,
                   scl = c(12, 12, 12),
                   locDiff = match(Sall[,1], sort(unique(Sall[,1]))),
                   from = 1, to = length(unique(Sall[,1])), 
                   pal = palUse16, 
                   meshColor = "faces", 
                   open3d = F)

# Module #17
colHighlight17 <- palPurple[round(Tree[which(Tree[,2] == 17), 3]*100, 0)]
palUse17 <- rep("#F0F0F0", length(unique(Sall[,1])))
palUse17[which(sort(unique(Sall[,1])) == 17)] <- colHighlight17
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,
                   scl = c(12, 12, 12),
                   locDiff = match(Sall[,1], sort(unique(Sall[,1]))),
                   from = 1, to = length(unique(Sall[,1])), 
                   pal = palUse17, 
                   meshColor = "faces", 
                   open3d = F)

# Module #18
colHighlight18 <- palGreen[round(Tree[which(Tree[,2] == 18), 3]*100, 0)]
palUse18 <- rep("#F0F0F0", length(unique(Sall[,1])))
palUse18[which(sort(unique(Sall[,1])) == 18)] <- colHighlight18
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,
                   scl = c(12, 12, 12),
                   locDiff = match(Sall[,1], sort(unique(Sall[,1]))),
                   from = 1, to = length(unique(Sall[,1])), 
                   pal = palUse18, 
                   meshColor = "faces", 
                   open3d = F)

# Module #19
colHighlight19 <- palGreen[round(Tree[which(Tree[,2] == 19), 3]*100, 0)]
palUse19 <- rep("#F0F0F0", length(unique(Sall[,1])))
palUse19[which(sort(unique(Sall[,1])) == 19)] <- colHighlight19
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,
                   scl = c(12, 12, 12),
                   locDiff = match(Sall[,1], sort(unique(Sall[,1]))),
                   from = 1, to = length(unique(Sall[,1])), 
                   pal = palUse19, 
                   meshColor = "faces", 
                   open3d = F)

# Module #20
colHighlight20 <- palRed[round(Tree[which(Tree[,2] == 20), 3]*100, 0)]
palUse20 <- rep("#F0F0F0", length(unique(Sall[,1])))
palUse20[which(sort(unique(Sall[,1])) == 20)] <- colHighlight20
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,
                   scl = c(12, 12, 12),
                   locDiff = match(Sall[,1], sort(unique(Sall[,1]))),
                   from = 1, to = length(unique(Sall[,1])), 
                   pal = palUse20, 
                   meshColor = "faces", 
                   open3d = F)

# Module #22
colHighlight22 <- palRed[round(Tree[which(Tree[,2] == 22), 3]*100, 0)]
palUse22 <- rep("#F0F0F0", length(unique(Sall[,1])))
palUse22[which(sort(unique(Sall[,1])) == 22)] <- colHighlight22
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,
                   scl = c(12, 12, 12),
                   locDiff = match(Sall[,1], sort(unique(Sall[,1]))),
                   from = 1, to = length(unique(Sall[,1])), 
                   pal = palUse22, 
                   meshColor = "faces", 
                   open3d = F)

# Module #23
colHighlight23 <- palRed[round(Tree[which(Tree[,2] == 23), 3]*100, 0)]
palUse23 <- rep("#F0F0F0", length(unique(Sall[,1])))
palUse23[which(sort(unique(Sall[,1])) == 23)] <- colHighlight23
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,
                   scl = c(12, 12, 12),
                   locDiff = match(Sall[,1], sort(unique(Sall[,1]))),
                   from = 1, to = length(unique(Sall[,1])), 
                   pal = palUse23, 
                   meshColor = "faces", 
                   open3d = F)
# Outermost circle completed.

# Module #21
colHighlight21 <- palRed[round(Tree[which(Tree[,2] == 21), 3]*100, 0)]
palUse21 <- rep("#F0F0F0", length(unique(Sall[,2])))
palUse21[which(sort(unique(Sall[,2])) == 21)] <- colHighlight21
newMesh21 <- addMesh(mMesh, 
                   x_center = (vert$rend[7]) * cos(vert$thetaend[7]),
                   y_center = (vert$rend[7]) * sin(vert$thetaend[7]),
                   z_center = 0,
                   scl = c(10, 10, 10),
                   locDiff = match(Sall[,2], sort(unique(Sall[,2]))),
                   from = 1, to = length(unique(Sall[,2])), 
                   pal = palUse21, 
                   meshColor = "faces", 
                   open3d = F)

# Module #13
colHighlight13 <- palRed[round(Tree[which(Tree[,2] == 13), 3]*100, 0)]
palUse13 <- rep("#F0F0F0", length(unique(Sall[,3])))
palUse13[which(sort(unique(Sall[,3])) == 13)] <- colHighlight13
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,
                   scl = c(10, 10, 10),
                   locDiff = match(Sall[,3], sort(unique(Sall[,3]))),
                   from = 1, to = length(unique(Sall[,3])), 
                   pal = palUse13, 
                   meshColor = "faces", 
                   open3d = F)

# Module #5
colHighlight5 <- palRed[round(Tree[which(Tree[,2] == 5), 3]*100, 0)]
palUse5 <- rep("#F0F0F0", length(unique(Sall[,6])))
palUse5[which(sort(unique(Sall[,6])) == 5)] <- colHighlight5
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,
                   scl = c(10, 10, 10),
                   locDiff = match(Sall[,6], sort(unique(Sall[,6]))),
                   from = 1, to = length(unique(Sall[,6])), 
                   pal = palUse5, 
                   meshColor = "faces", 
                   open3d = F)

# Module #2
colHighlight2 <- palRed[round(Tree[which(Tree[,2] == 2), 3]*100, 0)]
palUse2 <- rep("#F0F0F0", length(unique(Sall[,9])))
palUse2[which(sort(unique(Sall[,9])) == 2)] <- colHighlight2
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,
                   scl = c(10, 10, 10),
                   locDiff = match(Sall[,9], sort(unique(Sall[,9]))),
                   from = 1, to = length(unique(Sall[,9])), 
                   pal = palUse2, 
                   meshColor = "faces", 
                   open3d = F)

# Module #7
colHighlight7 <- palPurple[round(Tree[which(Tree[,2] == 7), 3]*100, 0)]
palUse7 <- rep("#F0F0F0", length(unique(Sall[,7])))
palUse7[which(sort(unique(Sall[,7])) == 7)] <- colHighlight7
newMesh7 <- addMesh(mMesh, 
                   x_center = (vert$rend[12]) * cos(vert$thetaend[12]),
                   y_center = (vert$rend[12]) * sin(vert$thetaend[12]),
                   z_center = 0,
                   scl = c(10, 10, 10),
                   locDiff = match(Sall[,7], sort(unique(Sall[,7]))),
                   from = 1, to = length(unique(Sall[,7])), 
                   pal = palUse7, 
                   meshColor = "faces", 
                   open3d = F)

# Module #8
colHighlight8 <- palPurple[round(Tree[which(Tree[,2] == 8), 3]*100, 0)]
palUse8 <- rep("#F0F0F0", length(unique(Sall[,7])))
palUse8[which(sort(unique(Sall[,7])) == 8)] <- colHighlight8
newMesh8 <- addMesh(mMesh, 
                   x_center = (vert$rend[15]) * cos(vert$thetaend[15]),
                   y_center = (vert$rend[15]) * sin(vert$thetaend[15]),
                   z_center = 0,
                   scl = c(10, 10, 10),
                   locDiff = match(Sall[,7], sort(unique(Sall[,7]))),
                   from = 1, to = length(unique(Sall[,7])), 
                   pal = palUse8, 
                   meshColor = "faces", 
                   open3d = F)

# Module #3
colHighlight3 <- palPurple[round(Tree[which(Tree[,2] == 3), 3]*100, 0)]
palUse3 <- rep("#F0F0F0", length(unique(Sall[,8])))
palUse3[which(sort(unique(Sall[,8])) == 3)] <- colHighlight3
newMesh3 <- addMesh(mMesh, 
                   x_center = (vert$rend[11]) * cos(vert$thetaend[11]),
                   y_center = (vert$rend[11]) * sin(vert$thetaend[11]),
                   z_center = 0,
                   scl = c(10, 10, 10),
                   locDiff = match(Sall[,8], sort(unique(Sall[,8]))),
                   from = 1, to = length(unique(Sall[,8])), 
                   pal = palUse3, 
                   meshColor = "faces", 
                   open3d = F)

# Module #11
colHighlight11 <- palGreen[round(Tree[which(Tree[,2] == 11), 3]*100, 0)]
palUse11 <- rep("#F0F0F0", length(unique(Sall[,7])))
palUse11[which(sort(unique(Sall[,7])) == 11)] <- colHighlight11
newMesh11 <- addMesh(mMesh, 
                   x_center = (vert$rend[22]) * cos(vert$thetaend[22]),
                   y_center = (vert$rend[22]) * sin(vert$thetaend[22]),
                   z_center = 0,
                   scl = c(10, 10, 10),
                   locDiff = match(Sall[,7], sort(unique(Sall[,7]))),
                   from = 1, to = length(unique(Sall[,7])), 
                   pal = palUse11, 
                   meshColor = "faces", 
                   open3d = F)

# Module #4
colHighlight4 <- palGreen[round(Tree[which(Tree[,2] == 4), 3]*100, 0)]
palUse4 <- rep("#F0F0F0", length(unique(Sall[,10])))
palUse4[which(sort(unique(Sall[,10])) == 4)] <- colHighlight4
newMesh4 <- addMesh(mMesh, 
                   x_center = (vert$rend[18]) * cos(vert$thetaend[18]),
                   y_center = (vert$rend[18]) * sin(vert$thetaend[18]),
                   z_center = 0,
                   scl = c(10, 10, 10),
                   locDiff = match(Sall[,10], sort(unique(Sall[,10]))),
                   from = 1, to = length(unique(Sall[,10])), 
                   pal = palUse4, 
                   meshColor = "faces", 
                   open3d = F)

# Done!

这是我的问题：虽然我指定将面部网格旋转 30 度（在我的代码中使用旋转矩阵

rotMat_y

），但不同位置的面部仍然定位不同，通常接近但不完全是 30 度。 我猜原因是只有绘制在 rgl 窗口中心的人脸会显示精确的 30 度旋转，而其他位置的人脸看起来不会旋转 30 度，即使它们被指定为旋转 30 度。 我要解决的问题是我想让所有的面都精确旋转30度。 x-y-z 维度中面的方向应该看起来相同

我想出的一个可能的解决方案是在 rgl 窗口中创建单独的子场景，这样我不仅可以将每个面旋转到所需的角度，而且每个面看起来真的好像旋转了 30 度。 @user2554330 再次提供了在 another of my question 中的 rgl 窗口中创建子场景的方法。但是，对于我的情况，我不知道如何在 rgl 窗口中适当地指定子场景的位置。我对带有

vert$rend

和

vert$thetaend

值的 rgl 子场景窗口的大小感到困惑。

虽然我的代码看起来很复杂，但问题归结为在 3D 极地树状图指定的适当位置创建子场景。如果创建单独的子场景的想法可以解决我的问题，关键应该是

x_center = (vert$rend[xx]) * cos(vert$thetaend[xx]) and y_center = (vert$rend[xx]) * sin(vert$thetaend[xx]),

的更改和

newSubscene3d(newviewport = c(, , , )

的正确规范。

除了创建子场景，也欢迎任何其他解决方案。