rggplot2geometryconvex-hull

How do I find the surface between lines based on intercept and slope which includes the origin?


I'm looking for a way to visualize the surface between a number of straight lines, which are defined in a dataframe through their intercepts and slopes. The surface I am looking for is the one that encloses the origin (0, 0).

The number of lines may vary (even though in the following simplified example I only have 6), and some of them may be redundant (i.e. they do not enclose the surface I am looking for because other lines are more constraining).

Let's take this simple dataframe:

df <- data.frame("Line" = c("A", "B", "C", "D", "E", "F"),
                 "Intercept" = c(4, 3, -2.5, -1.5, -5, -.5),
                 "Slope" = c(-1, 1, 2.4, -.6, -.8, .6))

Plotting these lines with ggplot2:

ggplot(data = df) +
  geom_vline(xintercept = 0) +
  geom_hline(yintercept = 0) +
  geom_abline(mapping = aes(intercept = Intercept, slope = Slope),
              colour = "red") +
  coord_cartesian(xlim = c(-6, 6), ylim = c(-6, 6))

Gives me the following output:

ggplot output for diagonal lines with geom_abline()

Basically I want to find intersections between the lines that enclose the origin (0, 0), disregarding the redundant one (the bottom left in this case, with intercept = -5 and slope = -0.8). Those 5 intersection points would then be used to plot the convex hull.

My main problem lies in finding the intersection points of the constraining lines (green points below) in order to be able to find the blue surface.

ggplot output for diagonal lines with geom_abline including intersection points and surface

QUESTION: Any suggestions on how to deal with this in R, ideally in a way that can be extended to larger dataframes (including more constraining and redundant lines)?

ADDITIONAL QUESTION: geom_abline() does not have a group aesthetic similar to geom_line(), which could be used to identify the line. Does anyone know a workaround to draw straight lines in ggplot2 based on slopes and intercepts (or two user-defined points of the line)?

Thanks in advance for any suggestions or (parts of) potential solutions!

UPDATE: In order to center the polygon around point (a,b) instead of the origin (0, 0), I have amended the original code (in particular the ìnnermost()`-function from @AllanCameron as follows:

innermost <- function(slopes, intercepts, a, b) {

  meetings <- function(slopes, intercepts) {
    meets_at <- function(i1, s1, i2, s2) {
      ifelse(s1 - s2 == 0, NA, (i2 - i1)/(s1 - s2))
    }
    xvals <- outer(seq_along(slopes), seq_along(slopes), function(i, j) {
      meets_at(intercepts[i], slopes[i], intercepts[j], slopes[j])
    })

    yvals <- outer(seq_along(slopes), seq_along(slopes), function(i, j) {
      intercepts + slopes *
        meets_at(intercepts[i], slopes[i], intercepts[j], slopes[j])
    })

    cbind(x = xvals[lower.tri(xvals)], y = yvals[lower.tri(yvals)])

  }

  xy <- meetings(slopes, intercepts)
  xy[,1] <- xy[,1] - a
  xy[,2] <- xy[,2] - b

  is_cut <- function(x, y, slopes, intercepts, a, b) {
    d <- sqrt(x^2 + y^2)
    slope <- y / x
    xvals <- (intercepts + slopes*a - b) / (slope - slopes)
    yvals <- xvals * slopes + (intercepts + slopes*a - b)
    ds <- sqrt(xvals^2 + yvals^2)
    any(d - ds > 1e-6 & sign(xvals) == sign(x) & sign(yvals) == sign(y))
  }

  xy <- xy[sapply(seq(nrow(xy)), function(i) {
    !is_cut(xy[i, 1], xy[i, 2], slopes, intercepts, a, b)
  }),]

  xy <- xy[order(atan2(xy[,2], xy[,1])),]
  
  xy[,1] <- xy[,1] + a
  xy[,2] <- xy[,2] + b
  
  as.data.frame(rbind(xy, xy[1,]))
}

Solution

  • Here's a solution that only requires a bit of geometry and algebra, using only base R. We can define a function, innermost that finds the x, y co-ordinates of the inner polygon and returns them in counter-clockwise order as a data frame. This allows you to create your ggplot by doing:

    ggplot(data = df) +
      geom_vline(xintercept = 0) +
      geom_hline(yintercept = 0) +
      geom_abline(mapping = aes(intercept = Intercept, slope = Slope),
                  colour = "red") +
      geom_polygon(data = innermost(df$Slope, df$Intercept), aes(x, y),
                   fill = "#99d9ea") +
      geom_point(data = innermost(df$Slope, df$Intercept), aes(x, y), 
                 color = 'green3') +
      coord_cartesian(xlim = c(-6, 6), ylim = c(-6, 6))
    

    enter image description here

    The function innermost is defined as follows:

    innermost <- function(slopes, intercepts) {
      
      meetings <- function(slopes, intercepts) {
        meets_at <- function(i1, s1, i2, s2) {
          ifelse(s1 - s2 == 0, NA, (i2 - i1)/(s1 - s2))
        }    
        xvals <- outer(seq_along(slopes), seq_along(slopes), function(i, j) {
          meets_at(intercepts[i], slopes[i], intercepts[j], slopes[j])
        })
        
        yvals <- outer(seq_along(slopes), seq_along(slopes), function(i, j) {
          intercepts + slopes *
            meets_at(intercepts[i], slopes[i], intercepts[j], slopes[j])
        })
        
        cbind(x = xvals[lower.tri(xvals)], y = yvals[lower.tri(yvals)])
        
      xy <- meetings(slopes, intercepts)
      is_cut <- function(x, y, slopes, intercepts) {
        d <- sqrt(x^2 + y^2)
        slope <- y / x
        xvals <- intercepts / (slope - slopes)
        yvals <- xvals * slopes + intercepts
        ds <- sqrt(xvals^2 + yvals^2)
        any(d - ds > 1e-6 & sign(xvals) == sign(x) & sign(yvals) == sign(y))
      }
      xy <- xy[sapply(seq(nrow(xy)), function(i) {
        !is_cut(xy[i, 1], xy[i, 2], slopes, intercepts)
      }),]
      xy <- xy[order(atan2(xy[,2], xy[,1])),]
      as.data.frame(rbind(xy, xy[1,]))
    }
    

    Explanation

    Firstly, it's straightforward to get the intersection of two straight lines. The formula of a straight line is given by y = mx + c, where m is the slope and c is the intercept. So if we have two different straight lines given by the equations y = m1x + c1 and y = m2x + c2 then they must meet where m1x + c1 = m2x + c2

    Rearranging this we get

    m1x - m2x = c2 - c1

    or

    (m2 - m1)x = c2 - c1

    which means the lines meet where

    x = (c2 - c1)/(m1 - m2)

    That is, if we have intercept1 and slope1 for the first line and intercept2 and slope2 for the second, then we can find the x value of the meeting point with this simple function:

    meets_at <- function(intercept1, slope1, intercept2, slope2) {
      ifelse(slope1 - slope2 == 0, NA, (intercept2 - intercept1)/(slope1 - slope2))
    }
    

    Note that if the lines are parallel, i.e. slope1 - slope2 == 0, they will not have a unique meeting point and this function should return NA

    We can use this function across all pairs of lines to get all the intersections by using outer:

    meetings <- function(slopes, intercepts) {
      
      xvals <- outer(seq_along(slopes), seq_along(slopes), function(i, j) {
                meets_at(intercepts[i], slopes[i], intercepts[j], slopes[j])
               })
      
      yvals <- outer(seq_along(slopes), seq_along(slopes), function(i, j) {
                intercepts + slopes *
                meets_at(intercepts[i], slopes[i], intercepts[j], slopes[j])
              })
      
      cbind(x = xvals[lower.tri(xvals)], y = yvals[lower.tri(yvals)])
    }
    

    For example, if we plot these points we will see we get all the intersections of the lines plotted:

    plot(seq(-6, 6), seq(-6, 6), type = 'n')
    for(i in seq(nrow(df))) abline(a = df$Intercept[i], b = df$Slope[i])
    xy <- meetings(df$Slope, df$Intercept)
    points(xy[,1], xy[,2], col = 'red')
    

    enter image description here

    This still leaves us the problem of finding only the innermost points that form the outline of your desired polygon. To do this, note what happens when we draw a line from the origin (0, 0) to each of the intersection points in the image above:

    xy <- xy[abs(xy[,1]) < 6 & abs(xy[,2] < 6),] # Remove intersections outside plot
    
    for(i in seq(nrow(df))) {
      segments(x0 = 0, y0 = 0, x1 = xy[,1], y1 = xy[,2], col = '#0000ff20')
    }
    

    enter image description here

    Note that the blue lines going from the origin to the innermost vertices (which are the ones we want to keep) are not crossed by any other lines. In other words, the vertices we want to discard are those from which you cannot draw a straight line to the origin without it being intersected by another line.

    We can therefore work out if there are any lines crossing the segments joining the origin to the intersections, and keep only those that are not crossed.

    We also need to arrange the final points in circumferential order. This is done by calculating the arctangent of the angle between the x axis and a line drawn to the point, which is just atan2(y, x). This gives us a number between -pi and pi, with which the points can be ordered counter clockwise starting from the 9 o-clock position:

    innermost <- function(slopes, intercepts) {
      xy <- meetings(slopes, intercepts)
      is_cut <- function(x, y, slopes, intercepts) {
        d <- sqrt(x^2 + y^2)
        slope <- y / x
        xvals <- intercepts / (slope - slopes)
        yvals <- xvals * slopes + intercepts
        ds <- sqrt(xvals^2 + yvals^2)
        any(d - ds > 1e-6 & sign(xvals) == sign(x) & sign(yvals) == sign(y))
      }
      xy <- xy[sapply(seq(nrow(xy)), function(i) {
        !is_cut(xy[i, 1], xy[i, 2], slopes, intercepts)
      }),]
      xy <- xy[order(atan2(xy[,2], xy[,1])),]
      as.data.frame(rbind(xy, xy[1,]))
    }
    

    We can use the above function to find the innermost polygon created by a group of lines. Plotting in base R, we can do:

    plot(seq(-6, 6), seq(-6, 6), type = 'n')
    for(i in seq(nrow(df))) abline(a = df$Intercept[i], b = df$Slope[i])
    xy <- innermost(df$Slope, df$Intercept)
    points(xy$x, xy$y, col = 'red')
    polygon(xy$x, xy$y, col = 'gray')
    

    enter image description here