Hi again,
I've had a go at Prof Ripley's suggestion (Strauss process, code
below). It works great for my limited purpose (qualitative drawing,
really), I can just add a few mild concerns,
- ideally the hard core would not be a fixed number, but a
distribution of sizes (ellipses).
- I could not quite understand the link between the window size,
intensity, and number of elements returned. Is there a simple
relation I've missed?
- in general, I just have no clue how rStrauss works, I guess Prof
Ripley's first reference cited in ?Strauss would be useful for that
matter as I could not find anything with google but a C code looking
a bit alike some genetic algorithm (death and birth of objects).
Am I on the right track in thinking that the idea is not so
dissimilar with the following approach:
1) create a first random distribution (Poisson, apparently)
2) give birth and death to new objects
3) carry on with this evolution until the sample distribution
respects the given intensity and interaction parameters (e.g, for
hard core, no overlap is permitted within the radius)
This looks a bit like the idea of an "electrostatic problem",
1) start with an initial configuration
2) define the total potential energy as a cost function to be
minimized, sum of a boundary (repulsive term), and short distance
repulsion
3) optimize the positions with respect to this objective function
I tried unsuccessfully this second approach (see below), but I'm sort
of concerned about the handling of all these parameters by optim().
rStrauss works beautifully anyway !
Thanks again,
baptiste
##############
# Code
##############
library(spatstat) # Strauss process
## ellipse function obtained from R mailing list ##
## could use library(car) as an alternative
ellipse <- function(x,y,width,height,theta, npoints=100,plot=F)
{
# x = x coordinate of center
# y = y coordinate of center
# width = length of major axis
# height = length of minor axis
# theta = rotation
# npoints = number of points to send to polygon
# plot = if TRUE, add to current device
# = if FALSE, return list of components
a <- width/2
b <- height/2
theta<-theta*pi/180
xcoord <- seq(-a,a,length=npoints)
ycoord.neg <- sqrt(b^2*(1-(xcoord)^2/a^2))
ycoord.pos <- -sqrt(b^2*(1-(xcoord)^2/a^2))
xx <- c(xcoord,xcoord[npoints:1])
yy <- c(ycoord.neg,ycoord.pos)
x.theta <- xx*cos(2*pi-theta)+yy*sin(2*pi-theta)+x
y.theta <- yy*cos(2*pi-theta)-xx*sin(2*pi-theta)+y
if(plot)
invisible(polygon(x.theta,y.theta,col="black"))
else
invisible(list(coords=data.frame(x=x.theta,y=y.theta),
center=c(x,y),
theta=theta))
}
getEllipse<-function(dataf,plot=TRUE){
ellipse(dataf[1],dataf[2],dataf[3],dataf[4],dataf[5],plot=plot)
}
N <- 200 # initial number (loosely related to the actual number of
elements)
##############################################
## positions from Strauss hard core process ##
##############################################
positions <- rStrauss(beta = 0.1, gamma = 0, R = 2, W= square(N/2))
# beta: intensity
# gamma: interaction. 0 : hard core (no overlap), 1: complete
randomness
# R: radius of interaction (size of the core)
# W: window
N2<-length(positions$x)# number of elements
rand.angles <- rnorm(N2,sd=45,mean=45) # ellipse angles
rand.a <- rnorm(N2,sd=0.4,mean=1) # ellipse semi-axes
rand.b <- rnorm(N2,sd=0.4,mean=1) #
par(bty="n")
plot(cbind(positions$x,positions$y) ,cex = rand.a,axes=F,
xlab="",ylab="",t="n") # just circles: type="p"
ell.pos<-cbind(positions$x,positions$y,rand.a,rand.b,rand.angles)
apply(ell.pos,1,getEllipse) -> b.quiet # draw ellipses
##############################################
##############################################
##############################################
## optimizing an "electrostatic potential" problem (not working)
##############################################
N <-200
p0 <- matrix(rnorm(2*N),ncol=2) # random starting point
# boundary repulsion potential #
x <- seq(-1.2,1.2,l=N)
delta <- 0.1
# 1D example #
y <- -1*( 1 - exp(-abs(x)/delta) - exp(abs(x)/delta))
#plot(x,y,t="l")
# 2D example
xy <- expand.grid(x = x, y=x)
z.x <- -1*( 1 - exp(-abs(xy[,1])/delta) - exp(abs(xy[,1])/delta))
z.y <- -1*( 1 - exp(-abs(xy[,2])/delta) - exp(abs(xy[,2])/delta))
z <- matrix(z.x+z.y, ncol=length(x))
#image(x=x,y=x,z=z)
boundary <- function(xy = p0[1,]){
z.x <- -1*( 1 - exp(-abs(xy[1])/delta) - exp(abs(xy[1])/delta))
z.y <- -1*( 1 - exp(-abs(xy[2])/delta) - exp(abs(xy[2])/delta))
z.x + z.y
}
dist.2d <- function(x = c(1, 0 ),y = c(0,1) , w = 0.1){
(drop((x - y) %*% (x - y)) + w^2 )^(-3/2)
}
obj <- function(p = p0){
bound <- sum(apply(p,1,boundary)) # boundary term
rel.positions <- expand.grid(x=p,y=p) # all relative positions
bound + sum(apply(rel.positions,1,dist.2d)) # cost
}
optim(p0,obj) # fails miserably, but this does not sound right anyway
_____________________________
Baptiste Auguié
Physics Department
University of Exeter
Stocker Road,
Exeter, Devon,
EX4 4QL, UK
Phone: +44 1392 264187
http://newton.ex.ac.uk/research/emag
http://projects.ex.ac.uk/atto
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
