Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Filled in R gradient curve

Tags:

r

3d

I have a curve I am using R to make (see code below):

library(rgl)

y = seq(-5,25,by=0.01)
x = seq(5,20,by=0.02)

sd = 0.3*x
NAs <- rep(NA, length(x)*length(y))
z <- matrix(NAs, length(x), byrow = T)
for(i in seq(1,length(x))) {
    for(j in seq(1,length(y))) {
        val = dnorm(y[j],mean=7.5,sd=sd[i])     
        z[i,j] = val
        if(z[i,j] < 0.02) {
            z[i,j] = NA
        }
    }
}

col <- rainbow(length(x))[rank(x)]        

open3d()
persp3d(x,y,z,color=col,xlim=c(5,20),ylim=c(5,10),axes=F,box=F,xlab="exp",ylab="obs",zlab="p")

And here's what it makes: enter image description here

If you rotate it a bit, you'd be able to see that this is a hollow tube type figure.

enter image description here

But I'm trying to make it be filled in (with the color gradient) so that it's not hollow. Imagine taking a slice at any location, and you'd get a 2D plane, not a 2D curve, if that makes sense. How can I do this?

like image 306
CodeGuy Avatar asked Jan 28 '13 20:01

CodeGuy


2 Answers

To fill a gap (a 2-d shape) in 3-d you should not use lines, since they are 1-d objects. Fill the gap with triagles or quads (flat objects with four corners) instead.

library(rgl)

y <- seq(-5,25,by=0.1)
x <- seq(5,20,by=0.2)
z <- outer(.3*x, y, function(my.sd, my.y) dnorm(my.y, mean=7.5, sd=my.sd))
z[z < .02] <- NA

col <- rainbow(length(x))[rank(x)]        
xn <- length(x)
yn <- length(y)

open3d()
persp3d(x, y, z, color=col, xlim=c(5,20), ylim=c(5,10), axes=F, box=F,
        xlab="exp", ylab="obs", zlab="p")
rgl.quads(rep(x[xn], (yn-1)*4),
          sapply(2:yn, function(i) y[i-c(0:1,1:0)]),
          sapply(2:yn, function(i) c(z[xn,i-0:1], 0, 0)),
          color=col[xn])

enter image description here

The outer and sapply commands might be confusing if you are not that familiar to R, but think of them as vectorized for loops. The outer call does an outer join of the coordinates to create all of z in one go and the sapplys extract the coordinates of the quads. The reason for avoiding for loops in R (or any other high level non-compiled language) is that they are terribly slow and also make the code bulky.

like image 132
Backlin Avatar answered Oct 20 '22 00:10

Backlin


The best way to do this, after spending lots of time figuring out something more elegant, is to manually add lines to fill the gap:

yy = seq(-5, 25, by=0.01)
xx = rep(5,length(yy))
sds = 0.7 * xx
val = rep(NA, length(xx))
for(i in seq(1,length(val))) {
    val[i] = dnorm(yy[i],mean=rep(7.5,length(xx[i])),sd=sds[i])
    t = 0.06
    if(val[i] > 0.02) {
        #val[i] = t
        lines3d(c(xx[i],xx[i]),c(yy[i],yy[i]),c(0.02,val[i]),color="red")
    }
}
like image 35
CodeGuy Avatar answered Oct 20 '22 00:10

CodeGuy