ROC Surface plot

Question

How Can I draw ROC surface plot in R by coding as shown in the figure?

library(scatterplot3d)
a1=1
a2=1
b1=4
b2=2
d1=9
d2=4
n=100
x1 <-x2 <-seq(.001,.95, length =n)
x=expand.grid(x1, x2) 
f=function(x1,x2){
rocss=(pnorm(qnorm(1-x2,d1,d2),b1,b2)-pnorm(qnorm(x1,a1,a2),b1,b2))
return(rocss)
}
dens=matrix(f(x[,1],x[,2]), ncol = sqrt(length(x[,1])))
s3d <-scatterplot3d(x1, x2, seq(min(dens), max(dens), length = length(x1)), type = "n", grid = FALSE, angle = -45, zlab = expression(f(x[1], x[2])), xlab = expression(x[1]), ylab = expression(x[2]), main = "Bivariate distribution",zlim=c(0,1),xlim=c(0,1),ylim=c(0,1)) 
for(i in length(x1):1) s3d$points3d(rep(x1[i], length(x2)), x2, dens[i,], type = "l")
for(i in length(x2):1) s3d$points3d(x1, rep(x2[i], length(x1)), dens[,i], type = "l")