[R-sig-Geo] minimum area bounding rectangles
Barry Rowlingson
b.rowlingson at lancaster.ac.uk
Wed Mar 16 00:47:31 CET 2011
Shazam:
minbb <- function(pts){
ch = chull(pts)
pts=pts[ch,]
np = nrow(pts)
pts=rbind(pts,pts[1,])
minbba = Inf ; bbth = NA; rotmin = NA
for(i in 1:np){
th = pi-atan2(pts[i+1,2]-pts[i,2],pts[i+1,1]-pts[i,1])
prot = rotxy(pts,th)
bba = diff(range(prot[,1])) * diff(range(prot[,2]))
if(bba < minbba){
xyb=cbind(
c(min(prot[,1]),max(prot[,1]),max(prot[,1]),min(prot[,1])),
c(min(prot[,2]),min(prot[,2]),max(prot[,2]),max(prot[,2]))
)
xyb=rbind(xyb,xyb[1,])
xyb= rotxy(xyb,-th)
rotmin=prot
minbba = bba
bbth = th
}
}
return(list(minbba=minbba,theta=th,pts=rotmin,box=xyb))
}
rotxy = function(pts,angle){
co = cos(angle)
si = sin(angle)
cbind(co * pts[,1] - si * pts[,2], si * pts[,1] + co * pts[,2])
}
And usage:
> pts=cbind(runif(60),runif(60))
> bb=minbb(pts)
> plot(bb$box,type="l")
> points(pts)
> pts2=rotxy(pts,pi/1.6)
> bb2=minbb(pts2)
> plot(bb2$box,type="l")
> points(pts2)
> Thanks Barry for the great tips. The MATLAB code will be a huge help -
> should be relatively straightforward to implement in R.
Wrote this myself with a helpful rotate function nabbed from
spatstat. Could be optimised a bit, I reckon.
Barry
More information about the R-sig-Geo
mailing list