[R] plotting an isosurface on a 3d plot
Duncan Murdoch
murdoch@dunc@n @end|ng |rom gm@||@com
Mon Jun 17 21:07:11 CEST 2019
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 using 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 using gmail.com <mailto:murdoch.duncan using 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
>
More information about the R-help
mailing list