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:
If you rotate it a bit, you'd be able to see that this is a hollow tube type figure.
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?
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])
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 sapply
s 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.
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")
}
}
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With