Area: Group Project (Group 18)

Pick's Theorem

The area $A$ of a simple polygon in 2D with integer coordinates can be computed by Pick's theorem[16], which is a simple method and just need to know the number of integer coordinates points on the boundary, $B$; and number of integer coordinates points inside the polygon, $I$. Then $$A = I + \frac{1}{2}B -1.$$

Also, such polygon is known as lattice polygon, and such points are lattice points.

By Pick's theorem, the area of polygon in Figure 6.1 is $A = 20 + 6 - 1 = 25$.

Also, by Monte Carlo method:

x <- c(1,3,5,9,8,2)
y <- c(2,3,3,-1,-3,0)
area_mc_method_v4(x,y,typepoly = FALSE, plotpoly = FALSE, point_method = "wind", mc_method = "antithetic")

## [1] 25.0008

which shows that Pick's theorem return the right answer for the given convex polygon.

Pick's theorem also works for polygon with holes, as the area only depends on the number of lattice points.

A proof of Pick's theorem can be found here by Tom Davis.

It is also possible to imply Pick's theorem into a R function, as we have a function to test whether a point is inside the polygon or not, we just need to test all lattice points within the minimum rectangle that cover the polygon. But the point in polygon function might not work well for the points on the boundary, because in the Monte Carlo method we don't actually care about a few points on the boundary.

But still, we can first work out the point on the boundary by using the line equation method: check whether a point lies on the domain of the line equation that forms the edge, and get rid of them from the set of lattice points. Then we can test the rest points with the point in polygon function, and count the number that inside the polygon, also count the point on the boundary in previous stage.

Before we start, let's write a function that will get rid of the redundant point of a polygon. That is, if there is three or more points in a line, it will remove the point in-between.

start_from_minx_v2 <- function(x,y) {
  x1 <- start_from_minx(x,y)$x
  y1 <- start_from_minx(x,y)$y
  n <- length(x)
  index <- c()
  grad <- gradient(x1,y1)
  grad[n+1] <- grad[1]
  for (i in 1:n) {
    if (grad[i] == grad[i+1]) {
      index <- c(index,(i+1))
      repeats <- TRUE
    }
    else repeats <- FALSE
  }
  if (repeats) return(list(x = x1[-index],y = y1[-index])) else return(list(x = x1,y = y1))
}

If that is the case, then the gradient will be the same between those points, which allow us to allocate such points.

The only problem is, if a point is the vertex, then it lies on two edges, and we might have the repeat result. But that's also fine, just collect the index of the points on the boundary, and remove the duplicate. Infinity gradient is not even an issue: remember in Monte Carlo (Part 2), there is a function to find the line equation that in the form of $x = f(y)$. Another possible problem is the stability of the point in polygon function, but it works perfectly well with most of the polygons, so that should not be a problem in this case.

Below is the code in R:

picks_theorem <- function(x,y,gridplot = TRUE) {
  x1 <- start_from_minx_v2(x,y)$x # Against three or more points in one of more line(s)
  y1 <- start_from_minx_v2(x,y)$y
  if (type_of_polygon_v2(x1,y1) == "Self-intersecting") return("Error: Self-intersecting polygon")
  n <- length(x1)
  grad <- gradient(x1,y1);grad[n+1] <- grad[1]
  xi <- min(x1):max(x1);yi <- min(y1):max(y1);xj <- c();yj <- c()
  x2 <- x1;x2[n+1] <- x2[1];y2 <- y1;y2[n+1] <- y2[1]
  for (i in 1:length(xi)) {
    for (j in 1:length(yi)) {
      yj <- c(yj,yi[j]);xj <- c(xj,xi[i])
    }
  }
  m <- length(xj)
  boundi <- c()
  for (i in 1:m) {
    for (j in 1:n) {
      if (!is.infinite(grad[j])) {
        if (linen_equ(x1,y1,j)(xj[i]) == yj[i]) {
          if (xj[i] <= max(x2[j],x2[j+1]) && xj[i] >= min(x2[j],x2[j+1])) {
            boundi <- c(boundi,i)
          }
        }
      } else {
        if (linen_equy(x1,y1,j)(yj[i]) == xj[i]) {
          if (yj[i] <= max(y2[j],y2[j+1]) && yj[i] >= min(y2[j],y2[j+1])) {
            boundi <- c(boundi,i)
          }
        }
      }
    }
  }
  boundary <- boundi[!duplicated(boundi)]
  xnb <- xj[-boundary];ynb <- yj[-boundary]
  k <- length(xnb)
  inside <- c()
  if (k != 0) {
    for (a in 1:k) {
      if (point_inside_v2(x1,y1,xnb[a],ynb[a],"wind") == 1) {
        inside <- c(inside,a)
      }
    }
  }
  b <- length(boundary)
  i <- length(inside)
  area <- i + b/2 - 1
  if (gridplot) {
    plot(x,y,pch="")
    polygon(x,y)
    points(xj,yj)
    points(xj[boundary],yj[boundary],col="yellow",pch=19)
    if (i != 0) points(xnb[inside],ynb[inside],col="purple",pch=19)
  }
  return(area)
}

Now try it with the polygon at the beginning:

x <- c(1,3,5,9,8,2)
y <- c(2,3,3,-1,-3,0)
picks_theorem(x,y)

## [1] 25

It also works well with a simple rectangle:

x <- c(1,3,5,5,2,1,1)
y <- c(0,0,0,4,4,4,2)
picks_theorem(x,y)

## [1] 16

By the way, Monte Carlo method is working as usual:

area_mc_method_v4(x,y,typepoly = FALSE, plotpoly = FALSE)

## [1] 16

Sometimes complicate function might not be working as usual for those basic polygons, so it is worth to test on them.

Let's see one more example:

xe7 <- c(5,4,7,6,5,6,6,3,4,5,3,2,3,1,3,2)
ye7 <- c(4,2,3,-1,1,0,2,1,1,-1,0,1,2,1,3,4)
picks_theorem(xe7,ye7)

## [1] 12

You can test the function at here.

→Section: Areas and Determinants

© Group 18