[R] Plotting irregular grid as image or persp

[David Forrest]

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/