I am trying to find a method for finding the point of intersection of 2 hyperbolae in R. Single-branch hyperbloae can be described by the following equation:
or
where (xi, yi) and (xj, yj) are the coordinates of 2 foci (i and j), and r is the difference in distance between a given point on the hyperbola (x, y) and each of the foci.
It seems the best way to visualise a hyperbola using R is to visualise contour lines (when contour level = 0) of a 3D hyperboloid, determined using the above equation and implementing it into a function
f1 <- function(xi, yi, xj, yj, x, y, r){
sqrt((xj - x)^2 + (yj - y)^2) - sqrt((xi - x)^2 + (yi - y)^2) - r
}
As an example, we can build a grid and visualise the contour level 0 for both hyperbolae:
library(tidyverse)
library(sf)
# sample points
tribble(
~name, ~x, ~y,
'a', 25, 25,
'b', 75, 25,
'c', 50, 75,
) %>%
{. ->> sample_points}
# sample hyperbolae
expand.grid(
x = seq(0, 100, length = 100),
y = seq(0, 100, length = 100)
) %>%
as_tibble %>%
mutate(
z1 = f1(xi = 25, yi = 25, xj = 75, yj = 25, x, y, r = 5), # i = point a, j = point b
z2 = f1(xi = 25, yi = 25, xj = 50, yj = 75, x, y, r = 5) # i = point a, j = point c
) %>%
{. ->> hyp_data} %>%
ggplot()+
geom_contour(aes(x, y, z = z1), breaks = 0, colour = 1)+
geom_contour(aes(x, y, z = z2), breaks = 0, colour = 2)+
geom_point(data = sample_points, aes(x, y, color = name), size = 3)+
coord_sf()
Currently, I can find the point of intersection by extracting the line data from the two geom_contour lines, using ggplot_build(), converting the line coordinates to an sf LINESTRING, and use st_intersection to find the point of intersection:
# extract the data from the ggplot using ggplot_build()
ggplot_build(
hyp_data %>%
ggplot()+
geom_contour(aes(x, y, z = z1), breaks = 0, colour = 1)+
geom_contour(aes(x, y, z = z2), breaks = 0, colour = 2)+
# geom_point(data = sample_points, aes(x, y, color = name), size = 3)+ # we dont need the point data in ggplot_build()
coord_sf()
) %>%
.$data %>% # keep only the data component
map(. %>% select(x, y) %>% as.matrix %>% st_linestring) %>% # keep coordinates, turn into a linestring
{. ->> lines1}
# make the lines an sf object
tibble(a = 1:2) %>%
mutate(
geom = st_sfc(lines1),
) %>%
st_as_sf %>% # then make the whole thing an sf object
# use st_intersection to find the point of intersection
st_intersection %>%
# then keep only the 'point' (exclude the original 'lines')
filter(
st_is(geom, 'POINT')
) %>%
{. ->> intersection_point} %>%
ggplot()+
geom_sf(colour = 'red', size = 5)+
geom_contour(data = hyp_data, aes(x, y, z = z1), breaks = 0, colour = 1)+
geom_contour(data = hyp_data, aes(x, y, z = z2), breaks = 0, colour = 2)+
geom_point(data = sample_points, aes(x, y, colour = name), size = 3)
The limitation of this method is that it relies upon the ability of ggplot's geom_contour to visualise the contour lines. When the absolute value of r approaches the distance between two points (dist between a-b = 50), the hyperbola narrows and eventually becomes a straight line 'behind' one of the points. Here, geom_contour can't build contour lines, and therefore no line data is created, so the intersection can't be found. See how the following plot only has 5 lines when there should be 6; z6 doesn't plot and causes a warning message:
expand.grid(
x = seq(0, 100, length = 100),
y = seq(0, 100, length = 100)
) %>%
as_tibble %>%
mutate(
# z = f1(x, y, nr)
z1 = f1(xi = 25, yi = 25, xj = 75, yj = 25, x, y, r = 5),
z2 = f1(xi = 25, yi = 25, xj = 75, yj = 25, x, y, r = 15),
z3 = f1(xi = 25, yi = 25, xj = 75, yj = 25, x, y, r = 25),
z4 = f1(xi = 25, yi = 25, xj = 75, yj = 25, x, y, r = 35),
z5 = f1(xi = 25, yi = 25, xj = 75, yj = 25, x, y, r = 45),
z6 = f1(xi = 25, yi = 25, xj = 75, yj = 25, x, y, r = 50),
) %>%
ggplot()+
geom_contour(aes(x, y, z = z1), breaks = 0, colour = 1)+
geom_contour(aes(x, y, z = z2), breaks = 0, colour = 2)+
geom_contour(aes(x, y, z = z3), breaks = 0, colour = 3)+
geom_contour(aes(x, y, z = z4), breaks = 0, colour = 4)+
geom_contour(aes(x, y, z = z5), breaks = 0, colour = 5)+
geom_contour(aes(x, y, z = z6), breaks = 0, colour = 6)+
geom_point(data = sample_points, aes(x, y, color = name), size = 3)+
coord_sf()
# Warning messages:
# 1: stat_contour(): Zero contours were generated
# 2: In min(x) : no non-missing arguments to min; returning Inf
# 3: In max(x) : no non-missing arguments to max; returning -Inf
z6 in theory should be represented by the dashed line below:
I have developed methods for creating these straight lines when contour lines can't be generated, but I would like to be able to use a "mathematical/algebraical" approach to finding the point of intersection of the two lines, and not rely on the spatial/mapping/countour approach within R.
I have explored options using functions such as uniroot() and solve(), but have had limited success, possibly due to my limited understanding of the fundamental mathematics and/or the fact that the equation describing the hyperbola has both x and y terms?
My current idea involves writing a pair of equations, with an equal right hand side of 0 (apologies for incorrect maths language?):
(or a trio of equations):
whereby as above, (xi, yi), (xj, yj) and (xk, yk) are the coordinates of three foci (i, j & k), and r is the range difference. xi, yi, xj, yj, xk, yk can be substituted for the coordinates of points a, b, c in the above example
Then, I think it should be possible to solve the equations using a similar process to that described here https://www.wikihow.com/Algebraically-Find-the-Intersection-of-Two-Lines, but I don't know how this translates into R code and/or existing R functions.
The aim of this is to determine the location of a transmitter (point of intersection), based on time-difference-of-arrival (TDOA) of transmissions received by receiving antennae (foci) and hyperbolic positioning principles and trilateration/multilateration (similar to how GPS works).
I will be dealing with thousands to millions of pairs of hyperbolae, so ideally this process is relatively fast and manageable by a regular computer. Also, not essential but hopefully achievable is the ability to 'scale up' a process to handle more than a pair/trio of equal equations (max ~10), in instances where the same transmission is received by more than three stations (eg 'number of satellites' when using GPS).
I would greatly appreciate absolutely any ideas or help anyone can offer, as I have been thinking about and working on trying to figure this out for quite some time. Thank you.
