Dihedral angles
stla opened this issue · 12 comments
Hello,
I think it would be nice to have the dihedral angles of the edges. This is what I use (with my package cgalMeshes) to plot the "exterior" edges (those forming a non flat dihedral angle), like this for example:
library(rgl)
septuaginta <- mesh3d(
rbind(
c(0.866025403784439, 0.5, 0),
c(0.5, 0.866025403784439, 0),
c(0, 1, 0),
c(-0.5, 0.866025403784439, 0),
c(-0.866025403784439, 0.5, 0),
c(-1, 0, 0),
c(-0.866025403784439, -0.5, 0),
c(-0.5, -0.866025403784438, 0),
c(0, -1, 0),
c(0.5, -0.866025403784439, 0),
c(0.866025403784438, -0.5, 0),
c(1, 0, 0),
c(0.75, 0.5, -0.433012701892219),
c(0.433012701892219, 0.866025403784439, -0.25),
c(-0.433012701892219, 0.866025403784439, 0.25),
c(-0.75, 0.5, 0.433012701892219),
c(-0.866025403784439, 0, 0.5),
c(-0.75, -0.5, 0.433012701892219),
c(-0.43301270189222, -0.866025403784438, 0.25),
c(0.433012701892219, -0.866025403784439, -0.25),
c(0.75, -0.5, -0.433012701892219),
c(0.866025403784439, 0, -0.5),
c(0.43301270189222, 0.5, -0.75),
c(0.25, 0.866025403784439, -0.433012701892219),
c(-0.25, 0.866025403784439, 0.433012701892219),
c(-0.43301270189222, 0.5, 0.75),
c(-0.5, 0, 0.866025403784439),
c(-0.43301270189222, -0.5, 0.75),
c(-0.25, -0.866025403784438, 0.43301270189222),
c(0.25, -0.866025403784439, -0.433012701892219),
c(0.433012701892219, -0.5, -0.75),
c(0.5, 0, -0.866025403784439),
c(0, 0.5, -0.866025403784439),
c(0, 0.866025403784439, -0.5),
c(0, 0.866025403784439, 0.5),
c(0, 0.5, 0.866025403784439),
c(0, 0, 1),
c(0, -0.5, 0.866025403784439),
c(0, -0.866025403784438, 0.5),
c(0, -0.866025403784439, -0.5),
c(0, -0.5, -0.866025403784438),
c(0, 0, -1),
c(-0.433012701892219, 0.5, -0.75),
c(-0.25, 0.866025403784439, -0.43301270189222),
c(0.25, 0.866025403784439, 0.433012701892219),
c(0.433012701892219, 0.5, 0.75),
c(0.5, 0, 0.866025403784439),
c(0.433012701892219, -0.5, 0.75),
c(0.25, -0.866025403784438, 0.43301270189222),
c(-0.25, -0.866025403784439, -0.43301270189222),
c(-0.433012701892219, -0.5, -0.75),
c(-0.5, 0, -0.866025403784439),
c(-0.75, 0.5, -0.43301270189222),
c(-0.433012701892219, 0.866025403784439, -0.25),
c(0.433012701892219, 0.866025403784439, 0.25),
c(0.75, 0.5, 0.43301270189222),
c(0.866025403784438, 0, 0.5),
c(0.75, -0.5, 0.43301270189222),
c(0.43301270189222, -0.866025403784438, 0.25),
c(-0.433012701892219, -0.866025403784439, -0.25),
c(-0.75, -0.5, -0.433012701892219),
c(-0.866025403784438, 0, -0.5)
),
triangles = cbind(
c(2, 43, 53),
c(2, 13, 23),
c(2, 24, 34),
c(34, 44, 2),
c(1, 13, 2),
c(2, 53, 3),
c(9, 58, 8),
c(33, 43, 2),
c(2, 23, 33),
c(2, 44, 54),
c(54, 1, 2),
c(18, 7, 8),
c(8, 28, 18),
c(38, 28, 8),
c(14, 24, 2),
c(2, 3, 14),
c(19, 9, 8),
c(8, 58, 48),
c(48, 38, 8),
c(8, 7, 59),
c(59, 49, 8),
c(8, 49, 39),
c(29, 19, 8),
c(8, 39, 29)
) + 1,
quads = cbind(
c(4, 5, 16, 15),
c(52, 61, 5, 4),
c(55, 56, 11, 0),
c(0, 11, 21, 12),
c(42, 51, 61, 52),
c(15, 16, 26, 25),
c(25, 26, 36, 35),
c(36, 37, 47, 46),
c(45, 46, 56, 55),
c(35, 36, 46, 45),
c(52, 53, 43, 42),
c(46, 47, 57, 56),
c(10, 11, 56, 57),
c(24, 25, 35, 34),
c(34, 35, 45, 44),
c(12, 13, 1, 0),
c(53, 52, 4, 3),
c(20, 21, 11, 10),
c(30, 31, 21, 20),
c(9, 10, 57, 58),
c(12, 21, 31, 22),
c(22, 23, 13, 12),
c(44, 45, 55, 54),
c(0, 1, 54, 55),
c(26, 27, 37, 36),
c(27, 28, 38, 37),
c(14, 15, 25, 24),
c(3, 4, 15, 14),
c(40, 41, 31, 30),
c(50, 51, 41, 40),
c(61, 60, 6, 5),
c(60, 61, 51, 50),
c(19, 20, 10, 9),
c(32, 41, 51, 42),
c(22, 31, 41, 32),
c(42, 43, 33, 32),
c(32, 33, 23, 22),
c(47, 48, 58, 57),
c(37, 38, 48, 47),
c(5, 6, 17, 16),
c(16, 17, 27, 26),
c(6, 7, 18, 17),
c(17, 18, 28, 27),
c(60, 59, 7, 6),
c(59, 60, 50, 49),
c(49, 50, 40, 39),
c(29, 30, 20, 19),
c(39, 40, 30, 29)
) + 1
)
library(cgalMeshes)
mesh <- cgalMesh$new(septuaginta)
vertices <- mesh$getVertices()
edges <- mesh$getEdges()
open3d(windowRect = 50 + c(0, 0, 512, 512), zoom = 0.8)
shade3d(septuaginta, color = "navy")
plotEdges(
vertices, edges, color = "gold"
)Ok. My example above is not not good because all edges are non-flat. A good example must have faces made of multiple triangles.
Back from my travels.
There might be other packages that do that computation. rgl is mainly about display rather than geometric computation, so it might not fit here unless you are specifically interested in doing it for mesh3d objects.
For those, I'm not sure what the function interface would be. Would you want something like
dihedralAngle <- function(mesh, triangle1, triangle2)
where triangle1 and triangle2 would be vectors of face numbers, and the output would compute angles between corresponding elements? As you mentioned, angles between quads aren't always well defined, but it might also make sense to compute angles between segments and triangles, or between two segments. I'm not sure what the function header should look like to do that.
Not sure. In my package cgalMeshes I have a function getEdges which return all the edges of a triangle mesh with their corresponding dihedral angle. Then I select those whose dihedral angle is less than 179° to plot them. But i don't know how to get the edges and their two incident faces with rgl.
That's not something rgl ever needs to do for display, so there are no explicit functions for it. But you could write one that works on triangles in a mesh object. For example,
library(rgl)
meshEdges <- function(mesh) {
triangles <- mesh$it
nt <- NCOL(triangles)
if (nt) {
# Get all the edges of the triangles
edges <- cbind(face = rep(seq_len(nt), 3),
v1 = c(triangles[1,],
triangles[2,],
triangles[3,]),
v2 = c(triangles[2,],
triangles[3,],
triangles[1,]))
ne <- nrow(edges)
# Sort the vertex numbers into increasing order in each edge
edges[, 2:3] <- t(apply(edges[, 2:3], 1, sort))
# Sort the edges by vertex numbers so we can see dups
edges <- edges[order(edges[,2], edges[,3]),]
# Find which edges appear more than once
dup <- duplicated(edges[,2:3])
# Drop edges which have no dups
singleton <- !dup &
c(!dup[-1], TRUE)
edges <- edges[!singleton,]
dup <- dup[!singleton]
if (any(dup[which(dup) - 1])) {
warning("some edges appear more than twice; only the first pair reported")
multiple <- dup &
c(dup[-1], FALSE)
edges <- edges[!multiple]
dup <- dup[!multiple]
}
edges <- cbind(edges[!dup, 2:3],
face1 = edges[!dup, 1],
face2 = edges[dup, 1])
}
edges
}
meshEdges(icosahedron3d())
Maybe that's what you already have, or were you using a different representation of the meshes? In any case, I don't think I'll add that to rgl at this point in time.
Didn't examine your code yet. Should be tricky as usual.
No that's not what I do, I use cgalMeshes and this is easy with CGAL. Otherwise there's a function in Rvcg which returns the edges with their incident faces I think (not sure). But I'm afraid cgalMeshes has been archived on CRAN, there are some sad problems with the RcppCGAL package.
I carefully read your code and I didn't understand yet how it works. How do you characterize the "exterior" edges in this code? (the edges with an acute dihedral angle). I'm afraid I miss something obvious.
Honestly I think such a function would be nice in rgl. The wire3d function does not do what is generally desired to plot the edges of an icosahedron for example.
The code I posted returned all edges, not exterior ones. To get those you would need to filter the matrix my function returned according to a test for the edge being exterior. As I said, code for that is not in rgl.
I don't understand your comment about the icosahedron. Could you expand on that? What do you want it to do that it doesn't do?
Ah ok.
Regarding my comment: wire3d will plot the edges of all triangles, not only the edges of the polyhedron (in case if the polyhedron has non-triangular and non-quad faces). But my example of the icosahedron was stupid: its faces are triangular. I like to plot a polyhedron with its "exterior" edges as lines or tubes (with cylinder3d). Or even only the exterior edges, without the faces:
Okay, I get it. You're right, rgl won't easily plot edges of a dodecahedron. It works with triangles and quads, and there's no data structure to represent a pentagon other than by building it up from those. Recognizing that an edge returned by my meshEdges function was interior corresponds to determining if two lines intersect: the line through the edge and the line through the other two points in the two faces. Rounding error would mean you'd need some tolerance, but it's equivalent to a linear system with 3 equations and 2 unknowns, so an explicit test should be fairly easy to work out.
Here's a modification of the previous meshEdges function that identifies interior edges using the idea mentioned above.
library(rgl)
meshEdges <- function(mesh) {
triangles <- mesh$it
nt <- NCOL(triangles)
if (nt) {
# Get all the edges of the triangles
edges <- cbind(face = rep(seq_len(nt), 3),
v1 = c(triangles[1,],
triangles[2,],
triangles[3,]),
v2 = c(triangles[2,],
triangles[3,],
triangles[1,]))
ne <- nrow(edges)
# Sort the vertex numbers into increasing order in each edge
edges[, 2:3] <- t(apply(edges[, 2:3], 1, sort))
# Sort the edges by vertex numbers so we can see dups
edges <- edges[order(edges[,2], edges[,3]),]
# Find which edges appear more than once
dup <- duplicated(edges[,2:3])
# Drop edges which have no dups
singleton <- !dup &
c(!dup[-1], TRUE)
edges <- edges[!singleton,]
dup <- dup[!singleton]
if (any(dup[which(dup) - 1])) {
warning("some edges appear more than twice; only the first pair reported")
multiple <- dup &
c(dup[-1], FALSE)
edges <- edges[!multiple]
dup <- dup[!multiple]
}
edges <- cbind(edges[!dup, 2:3],
face1 = edges[!dup, 1],
face2 = edges[dup, 1])
}
# Determine if we have an interior edge by seeing if the
# line through the edge intersects the line through the other
# two vertices
interior <- apply(edges, 1, function(edge) {
verts <- edge[1:2]
faces <- edge[3:4]
otherverts <- setdiff(unique(triangles[, faces]), verts)
abcd <- asEuclidean2(mesh$vb[, c(verts, otherverts)])
# X matrix should have columns d-c and a-b, y vector is a-c
y <- abcd[,1] - abcd[,3]
x <- cbind(abcd[,4] - abcd[,3], abcd[,1] - abcd[,2])
resids <- lsfit(x, y, intercept = FALSE)$residuals
sqrt(sum(resids^2)) < 1.e-8*sqrt(sum(y^2))
})
data.frame(edges, interior)
}
dodec <- dodecahedron3d()
edges <- meshEdges(dodec)
edges <- edges[!edges$interior,]
outline <- mesh3d(vertices = dodec$vb, segments = rbind(edges$v1, edges$v2))
shade3d(outline, lit = FALSE, col = "blue")
Created on 2023-10-01 with reprex v2.0.2
Thanks! It also works fine for the septuaginta:
IMHO, "interior" and "exterior" are not good words. "Interior" ("inside") could make someone think of a mesh edge with a grave acute angle, inside the convex hull but not on the convex hull.
Maybe another word to denote a "true" edge instead of "exterior"? A "ridge"?