R help - Aug 2004 - Plotting irregular grid as image or persp

Hi,
   I have an array of 2d node locations and an array triangles, and would
like to plot something like a image or persp.


An example of doing it with rgl is:

  library(ncdf)
  library(rgl)
  # wget http://www.maplepark.com/~drf5n/extras/teapot.nc
  teapot<-open.ncdf("teapot.nc")
  z<-get.var.ncdf(teapot,"tris")
  l<-get.var.ncdf(teapot,"locations")
  z<-as.vector(z)

 demo.teapot<-function(x=1){
  rgl.triangles(l[1,z],l[2,z],l[3,z])
  x<-readline("Hit any key to continue")
  rgl.pop()}

 demo.teapot() # produces a 3d rgl plot of the teapot

I wrote some rgl.* helper functions (
http://www.maplepark.com/~drf5n/cgi-bin/wiki.cgi?ElcircVisualization )
to drive rgl.points(), rgl.lines(), rgl.triangles(), rgl.quads() such that
I can call :

rgl.drf5n.polys(verts=l,conns=z) # plot points, wireframe, triangles,
# quads, depending on shape of conns)

but I'd really like to be able to plot an irregular grid with the
known triangles as specified colors, and I'm not sure
what would be the best package.

library(fields) extrapolates the image beyond the data
library(tripack) refits and plots its own data

What should I use? Suggestions?

Dave
-- 
 Dave Forrest
 drf at vims.edu                                    (804)684-7900w
 drf5n at maplepark.com                             (804)642-0662h
                                   http://maplepark.com/~drf5n/

Thanks Deepanyan,


On Sat, 28 Aug 2004, Deepayan Sarkar wrote:
...
> Yes, I think rgl would be the right tool for this. Even apart from the 3d
> acceleration issues, one of the problems with getting this in R would be
that R
> doesn't do raster graphics, and I don't think hidden surface
algorithms are
> very easy to implement in the R model.
...
> >
> >   library(ncdf)
> > #  library(rgl)
> >   teapot<-open.ncdf("teapot.nc")
> >    # wget http://www.maplepark.com/~drf5n/extras/teapot.nc
> >   edges<-get.var.ncdf(teapot,"tris")
> >   vertices<-get.var.ncdf(teapot,"locations")
> >
> >   xy<-vertices[c(1,2),] # this would be cooler with ?persp's
trans3d
> >
> >   plot(1:2,xlim=range(unlist(xy[1,])),
ylim=range(xy[2,]),type='n')
> >   apply(edges,2,function(x){polygon(t(xy[,x]))})
>
> I was playing around with this yesterday and got something similar (but
more
> general). I didn't send it to you because I wasn't sure if
that's what you
> wanted. Of course, I'm more familiar with the lattice version of
trans3d, so it
> uses that. There are 2 versions, one using grid, one with base graphics.
>
> As I said, there are glitches due to faulty drawing order of the triangles.
> Shading is also possible (as implemented in wireframe), but those
calculations
> are done in C code, so it would take a bit longer to carry over.
>
>
> library(grid)
> library(lattice)
>
>
> plotMesh.grid <-
>     function(l, z, rot.mat, dist = 0.1)
>     ## rot.mat: 4x4 transformation matrix
>     ## dist: controls perspective, 0 = none
> {
>     x <- ltransform3dto3d(l[,z], rot.mat, dist = dist)
>     id <- seq(length = ncol(x) / 3)
>     ord <- order(x[3, id * 3] + x[3, id * 3 - 1] +
>                  x[3, id * 3 - 2])
>     grid.newpage()
>     xscale <- range(x[1,])
>     yscale <- range(x[2,])
>     md <- max(diff(xscale), diff(yscale))
>     pushViewport(viewport(w = 0.9 * diff(xscale) / md,
>                           h = 0.9 * diff(yscale) / md,
>                           xscale = xscale,
>                           yscale = yscale))
>     id <-
>         as.vector(outer(1:3, (id[ord]-1) * 3, "+"))
>     grid.polygon(x = x[1,id],
>                  y = x[2,id],
>                  default.units = "native",
>                  gp = gpar(fill = "gray"),
>                  id = rep(id[ord], each = 3))
> }
>
>
>
> plotMesh.base <-
>     function(l, z, rot.mat, dist = 0.1, subset = TRUE)
>     ## rot.mat: 4x4 transformation matrix
>     ## dist: controls perspective, 0 = none
> {
>     x <- ltransform3dto3d(l[,z], rot.mat, dist = dist)
>     id <- seq(length = ncol(x) / 3)
>     ord <- order(x[3, id * 3] + x[3, id * 3 - 1] +
>                  x[3, id * 3 - 2])
>     xscale <- range(x[1,])
>     yscale <- range(x[2,])
>     plot(xscale, yscale, type = "n")
>     x <- cbind(x, NA)
>     id <-
>         as.vector(rbind(outer(1:3, (id[ord]-1) * 3, "+"),
>                         ncol(x)))
>     polygon(x = x[1,id],
>             y = x[2,id],
>             col = "gray")
> }
>
>
> rot.mat <- ltransform3dMatrix(list(y = -30, x = 40))
> plotMesh.grid(l, z, rot.mat, dist = 0)
> plotMesh.base(l, z, rot.mat, dist = 0)
Those are nice -- I did want a varying color however, and needed to
separate the calls to polygon:

 plotMesh.base<-function(vertices,edges,col,rot.mat=diag(4),dist=0.1,...){
  ## rot.mat a 4x4 homogeneous transformation matrix
  ## dist: controls perpective per lattice::ltransform3dto3d

  # rotate
  vertices<-ltransform3dto3d(vertices,rot.mat,dist)
    xscale <- range(vertices[1,])
    yscale <- range(vertices[2,])
    plot(xscale, yscale, type = "n")

  # find rough plot order

  ord<-order(apply(edges,2,function(x){sum(vertices[3,x])}))

  if (length(col) == 1){
    sapply(ord,function(x){
      polygon(vertices[1,edges[,x]],vertices[2,edges[,x]],col=col,...)})
  } else {
    sapply(ord,function(x){
      polygon(vertices[1,edges[,x]],vertices[2,edges[,x]],col=col[x],...)})
  }
  invisible(ord)
 }

rot.mat <- ltransform3dMatrix(list(z=45,y=30)) #;rot.mat

plotMesh.base(l,z,rot.mat=rot.mat,col=rainbow(dim(z)[2]),lty=0)


and the resultant image is:

   http://www.maplepark.com/~drf5n/images/teapot2.png


I posted some really rough notes at
http://www.maplepark.com/~drf5n/cgi-bin/wiki.cgi?RMeshVisualization

Dave
-- 
 Dave Forrest
 drf at vims.edu                                    (804)684-7900w
 drf5n at maplepark.com                             (804)642-0662h
                                   http://maplepark.com/~drf5n/