# NOT RUN {
## Comment: Simulate a Poisson sphere system,
## intersect, discretize and display results
library(unfoldr)
library(rgl)
library(plotrix)
spheres <- function(spheres, box=NULL, draw.box=FALSE, draw.axis=FALSE, ...) {
xyz <- apply(sapply(spheres, "[[", "center"),1,function(x) x)
sizes <- unlist(lapply(spheres,function(x) x$r))
spheres3d(xyz,radius=sizes,...)
if(draw.axis) {
axes3d(c('x','y','z'), pos=c(1,1,1), tick=FALSE)
#title3d('','','x','y','z')
}
axes3d(edges = "bbox",labels=TRUE,tick=FALSE,box=TRUE,nticks=0,
expand=1.0,xlen=0,xunit=0,ylen=0,yunit=0,zlen=0,zunit=0)
if(draw.box) {
x <- box$xrange[2]
y <- box$yrange[2]
z <- box$zrange[2]
c3d.origin <- rgl::translate3d(rgl::scale3d(rgl::cube3d(col="gray", alpha=0.1),
x/2,y/2,z/2),x/2,y/2,z/2)
rgl::shade3d(c3d.origin)
}
}
col <- c("#0000FF","#00FF00","#FF0000","#FF00FF","#FFFF00","#00FFFF")
#################################################################
## `beta` distribution for radii
#################################################################
lam <- 50
## parameter beta distribution (radii)
theta <- list("size"=list("shape1"=1,"shape2"=10))
# simulation bounding box
box <- list("xrange"=c(-1,4),"yrange"=c(-1.5,3.5),"zrange"=c(0,2))
## simulate and return full spheres system
## intersect with XZ plane and return full list of intersection profiles
## section plane xy
S <- simPoissonSystem(theta,lam,size="rbeta",box=box,type="spheres",
intersect="full",n=c(0,0,1),dz=0,pl=1)
# check resulting distribution
length(S$S)
summary(sapply(S$S,"[[","r"))
theta$size[[1]]/(theta$size[[1]]+theta$size[[2]]) # mean
## interior spheres:
## the ones which intersect one of the lateral planes (without top/bottom planes)
## showing spheres with color intersect
notIn <- sapply(S$S,function(x) !attr(x,"interior"))
spheres(S$S[notIn],box,FALSE,TRUE,color=col)
# not intersecting
In <- sapply(S$S,function(x) attr(x,"interior"))
spheres(S$S[In],box,color="gray")
## full sphere system
# open3d()
# spheres(S$S,box,FALSE,TRUE,color=col)
# planes3d(0,0,1,0,col="darkgray",alpha=1)
## draw intersections
sp <- S$sp
id <- sapply(sp,"[[","id")
open3d()
spheres(S$S[id],box,FALSE,TRUE,color=col)
planes3d(0,0,1,0,col="darkgray",alpha=1)
XYr <- t(sapply(sp,function(s) cbind(s$center[1],s$center[2],s$r)))
# centers
x <- XYr[,1]
y <- XYr[,2]
r <- XYr[,3]
win <- attr(sp,"win")
xlim <- win[[1]]
ylim <- win[[2]]
dev.new()
plot(x,y,type="n",xaxs="i", yaxs="i", xlab="x",ylab="y",xlim=xlim,ylim=ylim)
for(i in 1:nrow(XYr))
draw.circle(x[i],y[i],r[i],nv=100,border=NULL,col="black")
## digitize inersections
## `win` can also be omitted if not different
## from original itersection window (derived from the box)
W <- digitizeProfiles(sp, delta=0.01, win=win)
dim(W)
dev.new()
image(1:nrow(W),1:ncol(W),W,col=gray(1:0))
#################################################################
## Exact simulation of spheres with log normal radii
#################################################################
lam <- 100
## parameter rlnorm distribution (radii)
theta <- list("size"=list("meanlog"=-2.5,"sdlog"=0.5))
# simulation bounding box
box <- list("xrange"=c(0,5),"yrange"=c(0,5),"zrange"=c(0,5))
## simulate only 3D system
S <- simPoissonSystem(theta,lam,size="rlnorm",box=box,type="spheres",
intersect="original", perfect=TRUE, pl=101)
## show
open3d()
spheres(S[1:2000],box,TRUE,TRUE,color=col)
## check!
mean(log(sapply(S,"[[","r")))
sd(log(sapply(S,"[[","r")))
#################################################################
## Planar section
#################################################################
# planar section of exact simulated `rlnorm` sphere system:
# returns diameters for those section profiles
# which have their centers inside the intersection window
sp <- planarSection(S,d=2.5,intern=TRUE,pl=1)
hist(sp)
summary(sp)
mean(log(sp/2))
sd(log(sp/2))
#################################################################
## General intersection, all objects (inter=FALSE)
#################################################################
SP <- intersectSystem(S, 2.5, n=c(0,0,1), intern=FALSE, pl=1)
## show in 3D
id <- sapply(SP,"[[","id")
open3d()
spheres(S[id],box,TRUE,color=col)
planes3d(0,0,-1,2.5,col="black",alpha=1)
## 2D sections
XYr <- t(sapply(SP,function(s) cbind(s$center[1],s$center[2],s$r)))
# centers
x <- XYr[,1]
y <- XYr[,2]
r <- XYr[,3]
xlim <- box$xrange
ylim <- box$yrange
dev.new()
plot(x,y,type="n",xaxs="i", yaxs="i", xlab="x",ylab="y",xlim=xlim,ylim=ylim)
for(i in 1:nrow(XYr))
draw.circle(x[i],y[i],r[i],nv=100,border=NULL,col="black")
# digitize
W <- digitizeProfiles(SP, delta=0.01)
dev.new()
image(1:nrow(W),1:ncol(W),W,col=gray(1:0))
#################################################################
## Unfolding
#################################################################
ret <- unfold(sp,nclass=25)
## Point process intensity
cat("Intensities: ", sum(ret$N_V)/25, "vs.",lam,"\n")
## original diameters
r3d <- unlist(lapply(S,function(x) 2*x$r))
rest3d <- unlist(lapply(2:(length(ret$breaks)),
function(i) rep(ret$breaks[i],sum(ret$N_V[i-1]))))
op <- par(mfrow = c(1, 2))
hist(r3d[r3d<=max(ret$breaks)], breaks=ret$breaks, main="Radius 3D",
freq=FALSE, col="gray",xlab="r")
hist(rest3d, breaks=ret$breaks,main="Radius estimated",
freq=FALSE, col="gray", xlab="r")
par(op)
#################################################################
## Update intersection: find objects which intersect bounding box
#################################################################
idx <- updateIntersections(S)
sum(!idx) # objects intersecting
id <- which( idx != 1)
open3d()
spheres(S[id],box,TRUE,TRUE,color=col)
#################################################################
## user-defined simulation function
#################################################################
lam <- 50
theta <- list("p1"=-2.5,"p2"=0.5)
myfun <- function(p1,p2) { c("r"=rlnorm(1,p1,p2)) }
box <- list("xrange"=c(0,5),"yrange"=c(0,5),"zrange"=c(0,5))
S <- simPoissonSystem(theta,lam,rjoint=myfun,box=box,type="spheres",pl=1)
r <- sapply(S,"[[","r")
mean(log(r))
sd(log(r))
open3d()
spheres(S,box,TRUE,TRUE,color=col)
# }
Run the code above in your browser using DataLab