# the input region should be a SpatialPolygonsDataFrame class object
# this is the default when importing from GRASS / OGR

# min.rad.sep: number of minimum of radii (multiplier) between blocks
# min.rad.edge: number of minimum of radii (multiplier) between block center and edge of region

# side-effects:
# 1. plots the result
rnd_blocks <- function(region, n.blocks, block.edge.length, min.rad.sep=2, min.rad.edge=1)
{
# convert region to spatstat style 'owin'
reg.window <- as.owin.SpatialPolygons(SpatialPolygons(region@polygons))

# min distance between sampling blocks
# will also use this to enforce edge effect problems
# this is based on the target block size
# block.edge.length

# minimum circumscribing radius
min.r <- sqrt((block.edge.length/2)^2 + (block.edge.length/2)^2)


# compute region 'eroded by' 100m
reg.win.eroded <- erode.owin(w=reg.window, r=(min.r * min.rad.edge) )

# compute points > 100 meters apart
p.ssi <- rSSI(r=(min.r * min.rad.sep), n=n.blocks, win=reg.win.eroded)


# generate blocks centered on random points:
# TL -> TR -> BR -> BL -> TL
blocks.x <- cbind(
p.ssi$x - (block.edge.length/2), 
p.ssi$x + (block.edge.length/2), 
p.ssi$x + (block.edge.length/2), 
p.ssi$x - (block.edge.length/2), 
p.ssi$x - (block.edge.length/2)
)

blocks.y <- cbind(
p.ssi$y + (block.edge.length/2), 
p.ssi$y + (block.edge.length/2), 
p.ssi$y - (block.edge.length/2), 
p.ssi$y - (block.edge.length/2), 
p.ssi$y + (block.edge.length/2)
)



# convert to sp polygons:
polys <- apply(matrix(1:n.blocks), 1, function(i) {
p.coords <- data.frame(x=blocks.x[i, ], y=blocks.y[i, ])
Polygons(list(Polygon(p.coords)), ID=i)
}
)

# promote to 'SpatialPolygonsDataFrame' class:
# and copy the original spatial ref system data
SP <- SpatialPolygonsDataFrame(
SpatialPolygons(polys, proj4string=region@proj4string),
data=data.frame(block_num=1:n.blocks))


# some feedback:
plot(reg.window)		# spatstat class
points(p.ssi, pch=3)	# spatstat class
plot(SP, add=TRUE)		# sp class

# return results
SP
}