CodeGuy
CodeGuy

Reputation: 28905

Filled in R gradient curve

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?

Upvotes: 11

Views: 602

Answers (2)

Backlin
Backlin

Reputation: 14862

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.

Upvotes: 7

CodeGuy
CodeGuy

Reputation: 28905

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")
    }
}

Upvotes: 6

Related Questions