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:
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.
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,]))
}
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))
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 , where m is the slope and c is the intercept. So if we have two different straight lines given by the equations and then they must meet where
Rearranging this we get
or
which means the lines meet where
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')
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')
}
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')