[R-sig-Geo] minimum area bounding rectangles
Murray Richardson
murray_richardson at carleton.ca
Wed Mar 16 01:05:31 CET 2011
Wow, what can I say!? Thanks a million! I will try it out thoroughly
tomorrow.
Murray
On 15/03/2011 7:47 PM, Barry Rowlingson wrote:
> 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