plotting an isosurface on a 3d plot
See responses inline.
On 17/06/2019 1:18 p.m., ravi wrote:
Hi Duncan, Once again, I must thank you for your excellent support. I am getting to know things which I would have been very difficult by just reading the manual. I have followed your suggestion to use the qmesh3d command. I have uniquely defined the vertices. I am still stuck. I show below my code and the error message that I get. ##### isosurface on a 3d plot orig1 <- 0.3;radius1=0.2; theta <- seq(0,2*pi,pi/8) x1 <- orig1+radius1*cos(theta) y1 <- orig1+radius1*sin(theta) orig2 <- 0.5;radius2=0.2; theta <- seq(0,2*pi,pi/8) x2 <- orig2+radius2*cos(theta) y2 <- orig2+radius2*sin(theta) orig3 <- 0.3;radius3=0.1; theta <- seq(0,2*pi,pi/8) x3 <- orig3+radius3*cos(theta) y3 <- orig3+radius3*sin(theta) plot(x1,y1,type='b',col="red",xlim=c(0,1),ylim=c(0,1),asp=1) lines(x2,y2,type="b",col="blue") lines(x3,y3,type="b",col="green") z1 <- rep(0,length(x1)) z2 <- rep(0.5,length(x2)) z3 <- rep(1,length(x3)) ### use rgl package library(rgl) p1 <- cbind(x1,y1,z1) n <- length(x1) p1<- p1[-n,] plot3d(p1, col="red") p2 <- cbind(x2, y2, z2) p2<- p2[-n,] points3d(p2, col="blue") p3 <- cbind(x3, y3, z3) p3<- p3[-n,] points3d(p3, col="green") aspect3d(1,1,1) highlevel() vertices <- cbind(t(p1),t(p2),t(p3)) indices <- rbind(2:n, 1:(n-1),n + 1:(n-1), n + 2:n,2*n + 2:n, 2*n + 1:(n-1),n + 1:(n-1),n + 2:n)
That's too many rows. See the comment about indices below. It also includes values that are out of range. Your vertices matrix has 48 columns, but n is 17, so 2*n + 2:n goes up to 51.
cyl <- qmesh3d(vertices, indices, homogeneous = FALSE) str(cyl) shade3d(col = "grey",alpha=0.2) ## ########I get the following error message
> shade3d(col = "grey",alpha=0.2)
Error in UseMethod("shade3d") :
? no applicable method for 'shade3d' applied to an object of class
"character"
#######################
You forgot to put the mesh object in the call. It should be shade3d(cyl, col = "grey",alpha=0.2) or cyl %>% shade3d(col = "grey",alpha=0.2) (which is basically the same thing).
## the following works ind1 <- rbind(2,1,17,18,34,33,17,18) cyl1 <- qmesh3d(vertices, ind1, homogeneous = FALSE) str(cyl1) shade3d(addNormals(cyl1), col = "grey",alpha=0.2) ###################
> str(cyl)
List of 6 ?$ vb? ? ? ? ? ?: num [1:4, 1:48] 0.5 0.3 0 1 0.485 ... ?$ ib? ? ? ? ? ?: num [1:4, 1:32] 2 1 18 19 36 35 18 19 3 2 ...
> str(cyl1)
List of 6 ?$ vb? ? ? ? ? ?: num [1:4, 1:48] 0.5 0.3 0 1 0.485 ... ?$ ib? ? ? ? ? ?: num [1:4, 1:2] 2 1 17 18 34 33 17 18 Why does the former give an error message and not the latter?
addNormals(cyl1) gives a mesh object, so your second call had one.
Other questions: 1. Is there a better method of defining the indices? The qmesh3d requires that the length of indices be a sub-multiple of the number of rows. Can I get around this? For example, if I wanted to define the indices in the order c(2,1,17,33,34,18,2)? How flexible can I get in defining the indices vector?
The indices to qmesh3d() will be converted to a 4xn matrix, where each column gives the indices of the vertices forming one quadrilateral. You can pass a vector instead of a matrix and qmesh3d() will convert it to a 4 row matrix, but that's about all the flexibility you have.
2. The result with the shading command results in a surface with a white patch in the middle. Can this be avoided?
I'm not sure what you're referring to, but it might be the specular reflection from the surface. You can specify "lit=FALSE" to do away with all lighting, but that loses a lot of the 3D illusion. If you just want to get rid of the bright reflection, you can try shade3d(addNormals(cyl1), col = "grey",alpha=0.2, specular="black") to make it a non-reflective surface.
3. Just curious. When would one want to use homogenous coordinates in the definition of the vertices?
If they're more convenient :-). OpenGL does most calculations in homogeneous coordinates, because then affine transformations can be conveniently represented as matrices (i.e. linear transformations). There are also a few oddities: the sum in homogeneous coordinates is a weighted average in Euclidean coordinates.
4. If I want to paste separate images on the top and bottom planes, how would go about doing it? This is a case of texture mapping. I would really appreciate it if I get some tips.
It's easiest to plot the two surfaces separately. If you want to plot both at once, you need to glue the two images together into a single PNG, and use texture coordinates to choose which image is used for each vertex. It's probably a lot more trouble than just plotting the surfaces separately.
I am thankful for the rare previlige of directly getting help from the developer of the rgl package on this forum.
That's what this forum is for, and it's why most participants ask all responses to be sent to the list: if I answer your questions here, there's a permanent record which anyone can find via Google in the future. Duncan Murdoch
Thanking you, Ravi On Friday, 14 June 2019, 21:37:10 CEST, Duncan Murdoch <murdoch.duncan at gmail.com> wrote: On 14/06/2019 9:14 a.m., ravi wrote:
> Hi Duncan, > Thanks a lot for your help. You have really opened up a road for me to > to pursue further. > > Your later solution with cylinder3d is interesting. However, in my real > application, I have several parallel planes and very irregular contours > in each of them. So the quad3d command is more useful.? Are there more > efficient ways of using this without the for loops?
Yes, but you need to work out a longer vector of indices into the p1 and p2 vectors.? Depending on how you store your data that could be just tedious or really hard:? do you know how the vertices in p1 and p2 correspond?? If they are general curves with no connection, that can be a tricky question.
> > Other points where I would like to have help are : >? ? ? 1. I would like to smoothen the surface. But I have not been able > to figure out how I can use the addNormals command along with quad3d. > Are there other ways of improving the smoothness?
addNormals works on mesh objects, which can contain either quads or triangles.? If you've done the calculations to plot quads, you can create a qmesh3d object instead.? Create one matrix holding the coordinates of all vertices (ideally with no repetitions), and another matrix containing the indices of vertices which make up the quads. For example, starting with the code I posted before, this comes quite close: vertices <- cbind(t(p1), t(p2)) n <- nrow(p1) indices <- rbind(2:n, 1:(n-1), n + 1:(n-1), n + 2:n) cyl <- qmesh3d(vertices, indices, homogeneous = FALSE) shade3d(addNormals(cyl), col = "red") It's not quite right because vertices 1 and n have the same coordinates, as do vertices n+1 and 2*n; you really want to set it up so vertices are uniquely defined.
>? ? ? 2. As I said before, I can have 5 to 6 parallel planes. Sometimes, > the jump between two of them can be too high. So, I am thinking of using > an interpolation method at an intermediate plane. Do you have any > recommendations? I was thing of using a gam method where I make an > additive model for y ~ s(x) + s(z). Then predict y at height z given x. > I may have to figure out how I do this for closed curves (two or more > values of y for a given x). I just wonder if you can indicate generally > the route hat I should follow.
Nope, no idea on this one.? You might want to use the misc3d::contour3d function rather than working this all out for yourself, but it takes different inputs than you've described.
>? ? ? 3. How do I make the surface transparent? How do I add an alpha > value to the quad3d command?
That's a material property, so you can just put it in as part of the dots.? For the meshes, it would go into the shade3d() call in the same way, or the material argument to qmesh3d. Duncan Murdoch
> > Once again, I would like to thank you for your fantastic help. > Thanking you, > Ravi > > On Thursday, 13 June 2019, 22:52:07 CEST, Duncan Murdoch > <murdoch.duncan at gmail.com <mailto:murdoch.duncan at gmail.com>> wrote: > > > On 13/06/2019 4:32 p.m., Duncan Murdoch wrote: >? > On 13/06/2019 12:47 p.m., ravi via R-help wrote: >? >> Hi,I want to plot a surface joining a circle on a plane with another > circle (a little offset w.r.t. the first one) on a parallel plane. This > is a very simplified version of my problem. >? >> I will explain with the following code : >? >> # First, I plot the two circlesorig1 <- 0.4;radius1=0.3;theta <- > seq(0,2*pi,pi/8)x1 <- orig1+radius1*cos(theta)y1 <-
orig1+radius1*sin(theta)
>? >> orig2 <- 0.7;radius2=0.2;theta <- seq(0,2*pi,pi/8)x2 <- > orig2+radius2*cos(theta)y2 <- orig2+radius2*sin(theta) >? >> >
plot(x1,y1,type='b',col="red",xlim=c(0,1),ylim=c(0,1),asp=1)lines(x2,y2,type="b",col="blue")
>? >> #### the z coordinates on two parallel plnesz1 <- > rep(0,length(x1))z2 <- rep(1,length(x2)) >? >> >? >> I have tried to plot my 3d figure using the packages plot3D, misc3d, > rgl etc. All these packages require z as a matrix. In my case, all of > x,y and z are vectors. How can I plot this surface? In addition, I would > like to add similiar surfaces to the same plot. >? > >? > It's not true that rgl requires z as a matrix:? that's one
possibility,
>? > but there are others.? Here's one way to do what you want.
>? >
>? > # First, compute the two circles
>? > orig1 <- 0.4;radius1=0.3;theta <- seq(0,2*pi,pi/8)
>? > x1 <- orig1+radius1*cos(theta)
>? > y1 <- orig1+radius1*sin(theta)
>? > orig2 <- 0.7;radius2=0.2
>? > x2 <- orig2+radius2*cos(theta)
>? > y2 <- orig2+radius2*sin(theta)
>? >
>? > #### the z coordinates on two parallel plnes
>? > z1 <- rep(0,length(x1))
>? > z2 <- rep(1,length(x2))
>? >
>? > library(rgl)
>? > p1 <- cbind(x1,y1,z1)
>? > plot3d(p1, col="red")
>? > p2 <- cbind(x2, y2, z2)
>? > points3d(p2, col="blue")
>? >
>? > quads <- matrix(numeric(), 0, 3)
>? > for (i in seq_along(x1)[-1]) {
>? >? ? quads <- rbind(quads, p1[i-1,], p1[i,], p2[i,], p2[i-1,])
>? > }
>? > quads3d(quads, col = "green")
>? >
>? > There are more efficient ways to do this (the for loop would be
slow if
>? > you had bigger vectors), but this should give you the idea. > > Here's one of the more efficient ways: > > cylinder3d(rbind(c(0.4, 0.4, 0), c(0.7, 0.7, 1)), >? ? ? ? ? ? ? radius = c(0.3, 0.2), >? ? ? ? ? ? ? e1 = rbind(c(0,0,1), c(0,0,1)), sides=20) %>% >? ? ? ? ? ? ? addNormals() %>% >? ? ? ? ? ? ? shade3d(col = "green") > > > Duncan Murdoch