Hello,


Here is a code based on Greg's second idea (1st step: find the furthest
points of the polygon, 2nd step: find the circumscribed circle of a
triangle)

I wrote a function that computes the smallest circumscribing circle of a
triangle, and one that searches for the furthest vertices of a polygon.

The second function could help you find the median axis of a shape.

For the smallest circumscribing, I just took the formula from wikipedia.

You'll find a dummy example illustrating the search for the circle of a
triangle (2nd step).

You'll also find a spatial example including the search of the furthest two
points (1st step).

I've tested this code on a few countries from natural earth. It worked
quite well. Except sometimes, the polygon overlaps the circle.

Also, it performs the operation on the first entity. So, if the shape is of
type multipolygon, the algorithm will only work on the first polygon it
finds.

I think I now have a sufficient material to calculate the linearity index
but don't hesitate if you see something wrong in my code or have advice for
some algorithm improvements.


#####################
## DUMMY EXAMPLE ##
#####################

# http://en.wikipedia.org/wiki/Circumscribed_circle

findCircumTriangle <- function(tricoords) {
vert.x <- tricoords[,1]
vert.y <- tricoords[,2]

D <- 2*(vert.x[1]*(vert.y[2]-vert.y[3])
      + vert.x[2]*(vert.y[3]-vert.y[1])
      + vert.x[3]*(vert.y[1]-vert.y[2]))

circum.x <- ((vert.x[1]^2 + vert.y[1]^2)*(vert.y[2]-vert.y[3])
             +(vert.x[2]^2 + vert.y[2]^2)*(vert.y[3]-vert.y[1])
             +(vert.x[3]^2 + vert.y[3]^2)*(vert.y[1]-vert.y[2])) / D

circum.y <- ((vert.x[1]^2 + vert.y[1]^2)*(vert.x[3]-vert.x[2])
             +(vert.x[2]^2 + vert.y[2]^2)*(vert.x[1]-vert.x[3])
             +(vert.x[3]^2 + vert.y[3]^2)*(vert.x[2]-vert.x[1])) / D

dx <- (vert.x[1] - circum.x)
dy <- (vert.y[1] - circum.y)
r <- sqrt(dx^2 + dy^2)

return(list(x=circum.x, y=circum.y, radius=r))
}


coords <- cbind(c(1,0,3), c(.5,2,3))
circum <- findCircumTriangle(coords); circum

plot( 0:5, 0:5, type='n', xlab='', ylab='', asp=1)
polygon(coords)
theta <- seq(0,2*pi, length.out=250)
polygon( circum$x + circum$r*cos(theta),
         circum$y + circum$r*sin(theta), border='blue')


######################
## SPATIAL EXAMPLE ##
######################

library(rgdal)

countries <- readOGR("D:/DATAS/NATURAL
EARTH/ne_110m_admin_0_countries_lakes/ne_110m_admin_0_countries_lakes.shp",
"ne_110m_admin_0_countries_lakes")

country <- countries[24, ]
coords <- country@polygons[[1]]@Polygons[[1]]@coords

findMidBetweenFurthest <-  function(polycoords) {
  d <-  as.matrix(dist(polycoords))
  dmax <- max(d)

  first <- which(d==dmax, arr.ind=T)[1,1]
  second <- which(d==dmax, arr.ind=T)[1,2]
  mid <- apply(cbind(polycoords[first, ], polycoords[second, ]), 1,
FUN=mean)

  return(list(first=first , second=second, coords=mid, radius=abs(dmax/2)))

}

# 1st Step
mid <- findMidBetweenFurthest(coords)

# PLOTTING
theta <- seq(0,2*pi, length.out=250)
polygon( mid$coords[1] + mid$radius*cos(theta),
         mid$coords[2] + mid$radius*sin(theta), border='blue' )

d <- tail(as.matrix(dist(rbind(coords, mid$coords))), 1)
furthest <- which(d == max(d), arr.ind=T)[1,1]

# 2nd STEP
if (max(d) > mid$radius) {
  coordsTriangle <- rbind(coords[furthest, ], coords[mid$first, ],
coords[mid$second, ])
  circum <- findCircumTriangle(coordsTriangle)
} else {
  circum <- list(x=mid$coords[1], y=mid$coords[2], radius=mid$radius)
}

# PLOTTING
plot(country)
points(rbind(coords[mid$first, ], mid$coords, coords[mid$second, ]),
col="red")
theta <- seq(0,2*pi, length.out=250)
polygon( circum$x + circum$r*cos(theta),
         circum$y + circum$r*sin(theta), border='blue')


2013/1/26 Greg Snow <538280@gmail.com>

> The "dists" variable represents the distance between the points and the
> circle measured along the radius of the circle.  Positive distances mean
> that the point/vertex is inside the circle and negative distances mean that
> the point is outside of the circle.  Without some type of penalty the
> optimal circle would generally have some of the points inside and some
> outside (this would be a regression like model fitting the best circle
> instead of the best line).  Since we want all the points/vertices to be in
> the circle we add the penalty for any points outside the circle.  This is
> the penalty that worked for my examples, some other penalty may work
> better.
>
> For a more analytic solution see my second post.  That idea could possibly
> use some improvement as well, it is just what I thought of between waking
> up and actually getting out of bed (then refined somewhat while waiting for
> train into work to come).  It makes sense to me that it should work for
> most cases, maybe even all cases, but I have not worked out a formal proof
> of that.
>
>
> On Sat, Jan 26, 2013 at 3:57 AM, Berry Boessenkool <
> berryboessenkool@hotmail.com> wrote:
>
> >
> > Hi Greg,
> >
> > why did you add in particular
> > + 10*sum(dists<0)
> > in the function to be minimized?
> >
> > If I understand this correctly, this is the "punishment" of points
> outside
> > the current circle.
> > But if the points are very far apart (distances are large), then wouldn't
> > optim focus more on that part of the sum, instead of the number of points
> > outside as above?
> > Shouldn't the factor 10 be relative to the distance units?
> > Or is there a way to have optim only look at circles enclosing all
> vertices
> > of the convex hull?
> >
> > I was thinking there must be an analytical solution, but apparently not:
> > http://en.wikipedia.org/wiki/Smallest_circle_problem
> >
> > Berry
> >
> >
> >
> >
> > --
> > View this message in context:
> >
> http://r-sig-geo.2731867.n2.nabble.com/Finding-the-Circumscribing-Circle-of-a-polygon-tp7582372p7582388.html
> > Sent from the R-sig-geo mailing list archive at Nabble.com.
> >
> > _______________________________________________
> > 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
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]

