Thinking through this a little more I think I have a better algorithm to
start with:

Find the 2 vertices that are farthest from each other (the dist function
will compute all pairwise distances).  If you have highly concave polygons
or lots of vertices then using the chull (convex hull) may speed things up.

Find the circle that is centered on the point midway between those 2
vertices and has radius half the distance between them (it should intersect
the 2 points).  Check to see if all the other vertices are inside that
circle, if so then that should be the smallest enclosing circle.

If there are points outside the 1st circle then find the point that is
farthest outside the circle (the distance from the point to the center of
the circle is largest) and find the circle that goes through the 3 points
(this new one and the 2 used previously).  This new circle should include
all the vertices (and therefore the polygon) and I think it will be the
smallest possible.  If that circle is not the smallest enclosing circle
then at least it should be a good starting point for the optim routine.


On Fri, Jan 25, 2013 at 3:22 AM, Mathieu Rajerison <
mathieu.rajerison@gmail.com> wrote:

> Thanks a lot for your answer.
>
> I wasn't familiar with optimization methods.
>
> One question I asked myself was how to determine the optimal center
> coordinates and radius to start the optimization with.
>
> So, based on your code:
> - I used the centroid as the center
> - half the size of the biggest side of the bounding box as radius
>
> I got one problem. After optimizing, the circle still touched the shape.
> So, I made a loop that increased the radius until the shape didn't touch
> the circle.
>
> I simply translated your method to make it work on SpatialPolygons object.
>
> I used a sample from natural earth countries data as a sample to test your
> method and it worked quite well.
>
> Thanks again for your answer.
>
> The code, for the moment based on Greg Snow's solution (any advice
> appreciated)
>
> library(rgeos)
>
> myfun <- function(theta) {
>   x <- theta[1]
>   y <- theta[2]
>   r <- theta[3]
>
>   dists <- r - sqrt( (vert.x-x)^2 + (vert.y-y)^2 )
>   sum( dists^2 ) + 10*sum(dists<0)
> }
>
> spP = countries[13, ] # sampling a country from Natural Earth Data
> spP = SpatialPolygons(list(Srs1))
>
> # POLYGON COORDINATES
> vert.x = spP@polygons[[1]]@Polygons[[1]]@coords[, 1]
> vert.y = spP@polygons[[1]]@Polygons[[1]]@coords[, 2]
>
> # OPTIMIZATION
> coords <- coordinates(spP)
> ini.r <- min(diff(t(bbox(spP))))
> out <- optim(c(coords[1], coords[2], ini.r), myfun)
>
> # CIRCLE
> theta <- seq(0,2*pi, length.out=250)
> coordPol = cbind(out$par[1] + out$par[3]*cos(theta), out$par[2] +
> out$par[3]*sin(theta))
> coordPol = rbind(coordPol, coordPol[1, ])
> circle = SpatialPolygons(list(Polygons(list(Polygon(coordPol)), "pol")))
>
> # BUFFER CIRCLE UNTIL NO INTERSECTION
> i = 0
> inc = 0.01
> d = out$par[3] * inc
> while(!gContains(circle, spP)) {
>   d = d * (1 + inc)
>   circle <- gBuffer(circle, width = radius)
>   i = i+1
>   print(i)
> }
>
> # PLOTTING
> plot(spP)
> plot(circle, add=T)
>
> # LINEARITY INDEX
> lin = 1 - (gArea(spP)/gArea(circle))
>
>
> 2013/1/25 Greg Snow <538280@gmail.com>
>
>> Here is one possibility to get you close:
>>
>> vert.x <- c(0,0,1,1)
>> vert.y <- c(0,1,1,0)
>>
>> myfun <- function(theta) {
>> x <- theta[1]
>>  y <- theta[2]
>> r <- theta[3]
>>
>> dists <- r - sqrt( (vert.x-x)^2 + (vert.y-y)^2 )
>>  sum( dists^2 ) + 10*sum(dists<0)
>> }
>>
>> out <- optim( c(0,0,5), myfun )
>> out
>>
>>
>> plot( -1:2, -1:2, type='n', xlab='', ylab='', asp=1 )
>> polygon(vert.x, vert.y)
>>
>> theta <- seq(0,2*pi, length.out=250)
>> polygon( out$par[1] + out$par[3]*cos(theta),
>> out$par[2] + out$par[3]*sin(theta), border='blue' )
>>
>> tmp.x <- rnorm(25)
>> tmp.y <- rnorm(25)
>> tmp <- chull(tmp.x,tmp.y)
>>
>> vert.x <- tmp.x[tmp]
>> vert.y <- tmp.y[tmp]
>>
>> plot(vert.x, vert.y, xlim=c(-3,3), ylim=c(-3,3))
>> polygon(vert.x,vert.y)
>>
>> out <- optim( c(0,0,3), myfun )
>> out
>>
>> theta <- seq(0,2*pi, length.out=250)
>> polygon( out$par[1] + out$par[3]*cos(theta),
>> out$par[2] + out$par[3]*sin(theta), border='blue' )
>>
>>
>>
>> On Thu, Jan 24, 2013 at 7:19 AM, Mathieu Rajerison <
>> mathieu.rajerison@gmail.com> wrote:
>>
>>> I made a little mistake. The formula is 1 - area(shape) / area(circum),
>>> not 1 - area(circum) / area(shape)
>>>
>>> 2013/1/24 Mathieu Rajerison <mathieu.rajerison@gmail.com>
>>>
>>> > Hi,
>>> >
>>> >
>>> > I know spatstat::incircle to calculate the inscribed circle of a shape.
>>> >
>>> > I'd like to know if there was a function to calculate the
>>> circumscribing
>>> > circle of a polygon.
>>> >
>>> > The area of this circumscribing circle could then be used to calculate
>>> a
>>> > linearity index: 1 - area(circum) / area(shape)
>>> >
>>> > For a triangle, the circumcenter is at the intersection of the
>>> different
>>> > perpendicular bisectors of its sides and each corner would be at the
>>> same
>>> > distance to this center, distance being known as circumradius.
>>> >
>>> > For a polygon, it seems to me it's a lot more complicated.
>>> >
>>> > It's more a geometrical problem than a spatial one but if anyone has an
>>> > advice, I'd be glad of it..
>>> >
>>> >
>>> > Mathieu
>>> >
>>>
>>>         [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo@r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>>
>>
>> --
>> Gregory (Greg) L. Snow Ph.D.
>> 538280@gmail.com
>>
>
>


-- 
Gregory (Greg) L. Snow Ph.D.
538280@gmail.com

	[[alternative HTML version deleted]]

