Voronoi diagrams of moving points with Vorosketch

Vorosketch can draw Voronoi diagrams of points in the plane that have a speed (or more generally a magnitude of some sort) and a direction, applying one of various distance measures that take the speed and the direction of the sites into account. Note that these are not kinetic Voronoi diagrams; our goal here is not to keep a Voronoi diagram based on Euclidean distance up-to-date as the sites move. In contrast, the diagrams discussed on this page are static diagrams in which the definition of distance between a point and a site somehow involves the speed and the direction of the site.

In particular, the following distance measures are supported:

• catch@v, for any positive speed v
• push@a, for any acceleration a
• Dubins

All of these distance measures are invariant under translation and rotation of the coordinate system. Therefore, Vorosketch implements each distance measure d, where d is one the above, as a distance measure “oriented(@d)”, where @d is a distance measure that defines the distance to any point in the plane from a “site” that is assumed to be in the origin, directed upwards—leaving the speed or magnitude of the site as the only variable. The “oriented” function takes care of translating and rotating the coordinate system for each site such that these assumptions hold. On the rest of this page, we will therefore assume that a site is always located in the origin, directed upwards.

A Vorosketch input file with sites that have magnitude and direction starts with the number of sites, and then specifies each site by six numbers, separated by spaces:

2 x₁ y₁ x₂ y₂ w

The first number is always 2 (the number of defining points of the site). The location of the site is (x₁,y₁); its direction is given by the vector (x₂−x₁, y₂−y₁), normalised. its magnitude is the size of that vector, that is, ||(x₂−x₁, y₂−y₁)||. The last number is the weight of the site, and can simply be set to zero for an unweighted Voronoi diagram.

Here is an example input file, which we will take as input for the examples on this page:

5
2 -0.6 -0.6 -0.26 -0.6  0
2  0.6 -0.5  0.26 -0.5  0
2  0    0    0     0.14 0
2  0.2  0.2  0.5   0.5  0
2 -0.2  0.3 -0.5   0.6  0

(Implementing a more convenient input format is on the to-do list; I accept suggestions).

The @catch@v distance from a point p = (px, py) to a site in the origin, pointing upwards, is the minimum distance that an agent P, starting at p and moving with speed v, has to travel to meet an agent S that starts in the origin and moves upwards on the positive vertical axis at the site's specified speed w.

Let r be v/w. After time r·x/v has passed, S will be at coordinates (0, r·x·w/v) = (0, x), and could be met there by P if p lies on the circle with centre point (0, x) and radius r·x. If r > 1, then, as time passes, the circle grows in all directions and eventually reaches any point in the plane—in other words, P can eventually meet S, regardless of the starting point p. However, if r ≤ 1, then the circle's centre point moves upwards at least as fast as its radius grows, In that case, S can only be caught by P if P starts more or less in front of S.

To calculate the distance travelled by P, observe that both S and P will move to the meeting point m in a straight line. Thus, the origin o, p, and m form a triangle, in which ||m-p||, the distance travelled by P, equals r·||m-o|| = r·||m||. Let μ be the angle of this triangle at m. By the cosine rule we have the catch equation:

(r·||m||)² - ||m||² + 2·||m||·||p||·cos(μ) - ||p||² = 0,

where ||p||·cos(μ) = py. With the quadratic formula we find:

r·||m|| = r · ( -py ± √(py² + (r²-1)||p||²) ) / ( r²-1 ).

Indeed, if r < 1, a non-negative solution only exists if py ≥ 0 and py² + (r²-1)||p||² ≥ 0, that is, |px| / ||p|| ≤ r. In other words, p must lie in front of s, not too far to the side.

If r = 1, the quadratic formula does not give us a solution, but the catch equation simplifies to r·||m|| = ½·||p||²/py, which has a non-negative solution if and only if py > 0.

If r > 1, the quadratic formula always gives us exactly one positive solution, namely r·||m|| = r · ( -py + √(py² + (r²-1) ||p||² ) ) / (r²-1).

As an alternative to using the cosine rule, one can observe that m is the intersection of the positive vertical axis with the Apollonian circle of points whose distance to p is r times the distance to the origin. If r = 1, the Apollonian circle degenerates into the bisector of p and the origin. This results in the same formulas.

Unfortunately, these formulas are not numerically stable if r is very close to 1 (which could happen unintentionally due to rounding errors). Vorosketch currently solves this problem by using the formula for r = 1 whenever |r−1| ≤ 0.000001, in which case the Apollonian circle has radius over one million—enough so that for our purposes, it could be treated as a straight line. I am still looking for a more elegant solution that does not depend on a somewhat arbitrary threshold.

Here is a Voronoi diagram of the sites from the example input file given above:

vorosketch -m catch@0.36 -bczi 0.1 -w 1.1 -r 800 

Vorosketch also implements an unparameterised @catch distance. This is equivalent to the @catch@v distance when P and S always move at the same speed: the sites do not have specified speeds of their own. When using the unparameterised catch distance, the vectors specified in the input are only used to determine each site's direction; the length is ignored. The unparameterised catch distance also has a completely different interpretation.

The @push@a distance from a point p = (px, py) to a site (assumed to be in the origin, pointing upwards) is the minimum time in which an agent S, initially moving upwards from the site with the site's specified speed w, can reach an agent P at p, when S pushes itself into the right direction by applying a continuous acceleration of magnitude a.

Note that the time that S needs to reach P does not change if we add a vector -||p|| to each agent's location and a vector -w to each agent's speed, rotate the complete configuration by 180 degrees, and let ||P|| and ||S|| swap names. Thus, we can just as well assume that S starts in the origin and moves upwards with constant speed w, while P starts at ||p|| and accelerates from standstill to meet S. Obviously, in the transformed setting, P will move to the meeting point in a straight line, just like with the catch distance. The real difference with the catch distance is that P does not move with constant speed, but with constant acceleration (and is thus able to catch up with S eventually in all cases).

When S travels a distance ||m|| = w·t in time t = ||m||/w, the agent P now travels a distance ½·a·t². Thus, the cosine rule now gives us:

¼·a²·t⁴ - w²·t² + 2·py·w·t - ||p||² = 0,

and the distance between p and the site is defined as the smallest non-negative root of this quartic polynomial in t. Vorosketch implements a solution that might sometimes have issues with numerical stability, so be careful. The implementation also builds on unproven assumptions about which of the four roots of the quartic is needed in which cases. If you want to be safe, compile Vorosketch with -D CONSERVATIVE to insert a slightly slower snippet of code that does not build on these assumptions.

Here is a Voronoi diagram of the sites from the example input file given above:

vorosketch -m push@0.25 -bczi 0.1 -w 1.1 -r 800 

The @dubins distance from a site in the origin, pointing upwards, to a point p = (px, py), is the minimum distance that a “car” S, starting in the origin and pointing up, has to travel to reach p, when S can only move forward and can only change its direction by following curves whose curvature radius is at least the site's specified length r. I call this the Dubins distance because the path travelled is a Dubins path, generalised to end with any orientation. The structure of such paths has been analysed by Bui and Boissonat in their report "Accessibility region for a car that only moves forwards along optimal paths". My students Greta Günther and Felix Göhde showed me how to implement its computation.

If p lies inside the circle with centre point (-r,0) or (r,0) and radius r, then the path consists of a circular arc turning away from that circle, followed by a circular arc turning into the circle. Now let t be the centre of the other circle, let u be the centre of the appropriate tangent circle of the same radius that has p on the boundary (see Figure a) and let d be ||tp||/r. The length of the path travelled by S is now r·(α+β+2π−γ), where α = arctan(py/(|px|+r)); using the cosine rule on the triangle tpu, we find β = arccos(¼d + ¾ / d) and γ = arcsin(1¼ - ¼d²).

If p lies outside the circles with centre point (-r,0) and (r,0), then the path consists of a circular arc to turn into the right direction, followed by a straight line segment, see Figure b. The total length of the path is r·(α+β) + √(||tp||² − r²), where α = arctan((px−r)/py) and β = arccos(r/||tp||).

Correct implementation for all positions of p is achieved by strategic use of the atan2 function.

Note that points that are roughly in front of (above) the origin are easy to reach when thy lie outside the circles with centre point (-r,0) or (r,0) and radius r, and much harder when they lie just inside those circles. The car cannot make sharp enough turns to reach points inside those circles directly, and can only get there with a detour as in Figure a. As a result, the distance function is discontinuous on the front halves of those two circles. This can be seen clearly in the yellow region in the Voronoi diagram of the sites from the example input file given above:

vorosketch -m dubins -bczi 0.1 -w 1.5 -r 800