NOTE: This is a pretty old post, I wrote it but didn't publish it. Decided to publish it in the end.

I stumbled on this thing: The Reaktor Orbital Challenge and decided to give it a try (might not be online any longer by now).

If you haven't seen it before, or if you don't feel like heading over there to read it, then here's a recap.

Reaktor had a programming competition: you get a set of data representing satellites in orbit around Earth plus two sets of coordinates on the surface. The task is to route a fictional call from start to finish, through the satellites in orbit. You can win an Oculus Rift (pretty cool).

My first impression was that it looks a lot like a shortest-path-first problem, with one important caveat: you can only route between two satellites that have a clear line of sight. In other words, the Earth might get in between two potential hops and ruin it.

They also simplify things a bit by defining the Earth as a perfect sphere with the radius 6 371 km.

So, breaking it down I get a couple of partial problems:

  • From start, find the first satellite to route through.
  • From one satellite, find the next viable hop that will take me towards the finish.
    • Repeat until the next hop is the finish.

Breaking it down a little further, I get some more partial problems:

  • When looking for a potential next hop from any satellite, make sure there is a clear line of sight:
    • Drawing a line from one satellite to another, will it touch or intersect with the Earth?
    • Is there a maximum angle down (in the direction of the Earth) I can calculate to help limiting my choices?
  • All coordinates are given in latitude, longitude and altitude.
    • Would it be easier to switch to cartesian coordinates?
    • I think so. That would make things nice and clear (in my head).

TL;DR I've concocted a general solution based on some basic geometry and algebra

Breaking it down, this is my solution:

  • Convert all satellites and their terrestrial coordinates into points in 3D space with cartesian coordinates.
  • Use known math problems for finding out if a line intersects with a circle in 2D space to figure out if a vertex between two of the points have a clear line of sight.
  • Use known math problems for finding out if a point is above or below a plane in 3D space together with known math for finding the tangent plane of a sphere at a given point to figure out if a point in space is visible from a point on the surface of the sphere (Earth).
  • Use Dijkstra's algorithm for finding the shortest path to find the route from one terrestrial point, through a connected graph of nodes, ending at a final terrestrial point.

1.1. Switching coordinate space

Translating from latitude and longitude into simple () coordinates in 3D space can be done like this[1]:

Where is the radius of the Earth. Note that those will calculate points on the Earth, so we'll have to add the altitude we've got to move it out from the sphere.

Not to hard! In Ruby then:

def translate_coordinates(lat, lon, alt, earth_radius = 6371)
  nlat = lat.to_f * Math::PI / 180
  nlon = lon.to_f * Math::PI / 180
  r = alt.to_f + earth_radius
  [
    -(r) * Math.cos(nlat) * Math.cos(nlon),
    (r) * Math.sin(nlat),
    (r) * Math.cos(nlat) * Math.sin(nlon)
  ]
end

We'll keep that for later to use when importing the data snapshot.

1.2. Figuring out intersects with the Earth

I'm basing this on the same type of problem in 2D space[2], since that will translate straight up into 3D space[3] by simply adding a third spatial coordinate. Formulated in 2D space, this would be the problem of finding the intersection of a line between two points (satellites) and a circle (the Earth).

Here's a crappy drawing to illustrate the problem in 2D space:

A circle and two lines, one intersecting the circle.

We want to find the line from one point to another that doesn't intersect the sphere (or circle in the drawing above).

Given the Earth as a perfect sphere, we can use the cartesian equation to describe the sphere in 3D space.

Then, we can define our vector from to with all points along the line as

with being a real number .

Substituting for , and in the equation of the sphere we get the quadratic equation

which can be collected into

where

and then we can figure out the root(s) by calculating

I actually1 had to dig up my old TI-89 calculator for this, hadn't used that one in a looooong while.

In the final equation for calculating the root(s) , we can isolate the discriminant and since we have a quadratic equation that discriminant is in our case .

Then we can look at this discriminant to figure out if there are any intersect points between our line and sphere, or in someone else's words[4], replace with and ray with line to make it more applicable:

Depending on the values of , , and , the quadratic equation may have 0, 1, or 2 real solutions for . These three special cases correspond to three different scenarios where a ray may miss the sphere entirely (0 solutions), may graze the sphere tangentially at one point (1 solution), or may pierce the sphere, entering at one point and exiting at another (2 solutions). The value of the expression inside the quadratic solution's square root (called the radicand) determines which of the three special cases we counter:

  1. : no solutions
  2. : one solution
  3. : two solutions

Our function will regard the first case as two nodes having a clear line of sight or true and the other as there is a sphere in the way or false.

Now, if we translate all this into code, we'll get a function that returns a point of intersection. This can be used to make a function that simple replies truthy or falsey to the question can this point see this other point, or is there a planet in the way?

Ok, code then:

def line_of_sight?(p0, p1, from_surface = false, earth_radius = 6371)
  # Assuming p0 and p1 are objects with accessors :x, :y and :z

  # First calculate components of a quadratic equation derived from the math
  # we're dealing with.
  # The quadratic equation derived from [1] and [3] in terms of some
  # rational *t* is in our case "at^2 + bt + c = 0" with *a*, *b* and *c*
  # being the following:
  la = a(p0, p1)                # local calculated a
  lb = b(p0, p1)                # local calculated b
  lc = c(p0, p1, earth_radius)  # local calculated c

  clear_line = false  # assume there is no line of sight

  # We don't need to calculate the roots to see *if* there are any
  # intersect points. We can stop and look att the discriminant to determine
  # that. Since we have a quadratic equation in terms of the root *t*, the
  # discriminant in our case is "b^2 - 4ac". From [2] and [3] we get that
  # the following three cases exist based on the value of the discriminant:
  # b^2 - 4ac is positive (two intersections)
  if lb ** 2 - (4 * la * lc) > 0.0
    clear_line = false  # no line of sight
  end
  # b^2 - 4ac is zero (one intersect point, line is tangent)
  if lb ** 2 - (4 * la * lc) == 0.0
    clear_line = false  # no line of sight
  end
  # b^2 - 4ac is negative (no intersect)
  if lb ** 2 - (4 * la * lc) < 0.0
    clear_line = true   # clear line of sight between p0 and p1
  end
  clear_line
rescue
  # Something went amiss, just return false. If there were problems like
  # division by zero then we'll assume there weren't any intersection points.
  return false
end

# Some helpers below
private

def a(p0, p1)
  (p1.x-p0.x) ** 2 + (p1.y-p0.y) ** 2 + (p1.z-p0.z) ** 2
end

def b(p0, p1)
  2 * (((p1.x - p0.x) * p0.x) + ((p1.y-p0.y) * p0.y) + ((p1.z-p0.z) * p0.z))
end

def c(p0, p1, earth_radius = 6371)
  p0.x ** 2 + p0.y ** 2 + p0.z ** 2 - earth_radius ** 2
end

1.3. Finding nodes visible from the surface

The start and end points are a bit special since they are on the surface of the sphere so the algorithm above won't work that well. Instead, I'm going to look at it as finding the distance from a point on a plane to another point above that plane. The plane being the tangent plane of the sphere at the given point of its surface (the start or end point).

We can use the known position of the point on the surface of the sphere and the fact that our sphere is centered at to find the normal vector of the tangent plane as which gives us the tangent plane as [5].

From that, we can derive the distance from our given point[6] to some other point as:

If that distance is then the point is above the plane (in the direction of the normal vector), and if it's then it's below the plane (in the opposite direction of the normal vector).

If we convert this to Ruby code, we get something like this:

def find_nodes_above_tangent_plane_at_point(n0, nodes, earth_radius = 6371)
  above = []

  nodes.each do |n1|
    dist = (n0.x * n1.x + n0.y * n1.y + n0.z * n1.z - earth_radius ** 2) / Math.sqrt(n0.x ** 2 + n0.y ** 2 + n0.z ** 2)

    # Nodes "above" the tangent plane will have a positive distance to the
    # plane, above just means on the side towards which the normal vector
    # points.
    above << n1 if dist > 0 # n0 is above tangent plane
  end
  above   # return the list of node above the tangent plane
end

1.4. Finding a path

Now that we have a basic set of nodes in some space, each with coordinates, we can start thinking about this as a proper shortest-path-first problem. A good algorithm here would be something like Dijkstra's algorithm which pretty much does exactly what we want, if we add a constraint to stop the search once the target node is reached (the goal coordinates on the Earth). We can use the line of sight-thing to define the distance used in Dijkstra's algorithm to make the algorithm not pick those paths but that would probably result in those paths being picked anyway if no other exist (the case where a route does not actually exist), or we can use it when defining the graph of nodes (the satellites in orbit) so that certain nodes simply aren't connected to some other node (probably the best use of that data) and define all distances as the actual distances between two nodes (simple algebra when we have the cartesian coordinates to work with).

1.4.1. Distance

The distance between two points in 3D space is simply defined as:

or in Ruby:

module Helpers
  def node_dist(p0, p1)
    Math.sqrt((p1.x - p0.x) ** 2 + (p1.y - p0.y) ** 2 + (p1.z - p0.z) ** 2)
  end
end
1.4.2. Building a graph

Since I'm using Ruby, I'm going to set up some objects to represent the satellites/nodes in the graph. This way, I can let each node in the graph hold a list of reachable neighbours together with the distance to them. You could probably do this in any number of ways but this works for me.

So the satellites will look something like this:

class Node
  include Helpers

  attr_accessor :x, :y, :z
  attr_accessor :identifier

  def initialize(identifier, x, y, z)
    @identifier = identifier
    @x = x
    @y = y
    @z = z
  end

  def find_nodes_above_tangent_plane(nodes, earth_radius = 6371)
    vertexes = []
    nodes_above = find_nodes_above_tangent_plane_at_point(self, nodes, earth_radius)
    nodes_above.each do |n|
      vertexes << Vertex.new(self, n, distance(n))
    end
    vertexes
  end

  def setup_neighbours(nodes)
    vertexes = []
    nodes.each do |node|
      if in_line_of_sight?(node) && node.identifier != @identifier
        vertexes << Vertex.new(self, node, distance(node))
      end
    end
    vertexes
  end

  def distance(node2)
    node_dist(self, node2)
  end

  def in_line_of_sight?(node2)
    line_of_sight?(self, node2)
  end
end

Where the vertices looks something like this:

class Vertex
  include Helpers

  attr_accessor :endpoints, :distance

  def initialize(node1, node2, distance)
    @endpoints = [node1.identifier, node2.identifier]
    @distance = distance
  end

  # This needs to be available in order to run +#uniq+ on arrays with Vertex
  # objects.
  def eql?(other)
    endpoints_match = true
    @endpoints.each do |e|
      unless other.endpoints.member?(e)
        endpoints_match = false
      end
    end
    endpoints_match && @distance == other.distance
  end

  # This needs to be available in order to run +#uniq+ on arrays with Vertex
  # objects.
  def hash
    @endpoints.sort.hash
  end
end

Then the graph as something like:

class Graph
  attr_accessor :nodes, :vertexes

  def initialize
    @nodes = []
    @vertexes = []
  end

  def add_node(node)
    @nodes << node
  end

  def add_nodes(nodes)
    nodes.each do |node|
      add_node(node)
    end
  end

  def add_vertexes(vertexes)
    vertexes.each do |v|
      @vertexes << v
    end
    @vertexes.uniq!
  end

  def setup_neighbours
    @nodes.each do |node|
      neighbours = node.setup_neighbours(@nodes)
      neighbours.each do |n|
        @vertexes << n
      end
    end
    @vertexes.uniq!
  end
end
1.4.3. Traversing the graph

Going through the graph, we can use Dijkstra's algorithm with something like this in the Graph class (I'm not going to go through the details of Dijkstra's algorithm, there's a perfectly good Wikipedia page for that):

class Graph
  def traverse_graph(start_node, end_node)
    q = []
    dist = {}
    prev = {}

    @nodes.each do |node|
      dist[node.identifier] = Float::INFINITY # Unknown distance from source to v
      prev[node.identifier] = nil             # Previous node in optimal path from source
      q << node.identifier                    # All nodes initially in Q (unvisited nodes)
    end
    dist[start_node.identifier] = 0           # Distance from source to source

    while !q.empty?
      # Source node will be selected first
      nearest_node_identifier = dist.select{ |id,dist| q.member?(id) }.sort_by{ |identifier,distance| distance }.first.first
      # u ← vertex in Q with min dist[u]
      u = @nodes.select{ |node| node.identifier == nearest_node_identifier }.first

      if u.identifier == end_node.identifier
        s = []
        u = end_node
        while !prev[u.identifier].nil?
          s.unshift(u.identifier)
          u = @nodes.select{ |node| node.identifier == prev[u.identifier] }.first
        end
        s.unshift(u.identifier)
        return {
          seed: @seed,
          path: s
        }
      end
      # remove u from Q
      q.reject!{ |name| name == u.identifier }

      neighbour_vertexes = @vertexes.select{ |vertex|
        vertex.endpoints.member?(u.identifier)
      }

      neighbour_identifiers = neighbour_vertexes.map{ |v|
        v.endpoints
      }.flatten.reject{ |n|
        n == u.identifier || !q.member?(n)
      }.uniq

      neighbour_identifiers.each do |nv|
        v = @nodes.select{ |n| n.identifier == nv }.first
        alt = dist[u.identifier] + u.distance(v)
        if alt < dist[v.identifier]
          dist[v.identifier] = alt
          prev[v.identifier] = u.identifier
        end
      end
    end
    return [dist, prev]
  end
end

2. Wrapping up

Ok so I think we have everything we need to make this work, I'm thinking that the algorithm will go something like this:

def find_a_route(file="data/snapshot.txt")
  reader = SnapshotDataReader.new
  graph = reader.load(file)
  start_node = graph.nodes.select{ |n| n.identifier == "START" }.first
  end_node = graph.nodes.select{ |n| n.identifier == "END" }.first

  graph.traverse_graph(start_node, end_node)
end

solution = find_a_route("data/snapshot.txt")

So, this seems to be working just fine! Downloading a couple of data snapshots and verifying the calculated route tells me it's working as supposed:

Solved!

Fun challenge all in all! Guess all that's left now is hoping to win that Rift headset 😊 (Hint: I didn't, still fun though.)

3. References

  1. Converting latitude and longitude to x y z coordinates in Partiview. From COSMUS, University of Chicago Department of Astronomy and Astrophysics. http://astro.uchicago.edu/cosmus/tech/latlong.html
  2. Weisstein, Eric W. Circle-Line Intersection. From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Circle-LineIntersection.html
  3. Hinson, Anthony General Algorithms Section 2.1.4: Intersection of Parametric Line and Sphere. http://www.ahinson.com/algorithms_general/Sections/Geometry/IntersectionOfParametricLineAndSphere.pdf
  4. Cross, Don Fundamentals of Ray Tracing Chapter 6: Sphere Intersections. http://cosinekitty.com/raytrace/chapter06_sphere.html
  5. Nikunj (http://math.stackexchange.com/users/287774/nikunj), Tangent plane to sphere, URL (version: 2016-04-23): http://math.stackexchange.com/q/1517146
  6. Weisstein, Eric W. Point-Plane Distance. From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Point-PlaneDistance.html

  1. I can barely type the word actually without thinking about this one: https://www.facebook.com/buzzfeedadam/posts/1151999248184613 by the generally awesome Adam Ellis