[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