vorosketch

Vorosketch is a simple tool that sketches small Voronoi diagrams as a bitmap in a few seconds. Vorosketch's raison d'être is that it is easy to adapt to different distance measures, regardless of whether we understand the geometry of the bisectors yet. Thus we can use Vorosketch for explorations and illustrations; it supports several novel distance measures. Vorosketch draws:

- additively/substractively or multiplicatively/divisively weighted
- first-order, second-order, second-closest-site and farthest-site Voronoi diagrams under the
*Euclidean*(L2),*squared Euclidean*(power diagrams),*Manhattan*(L1) metrics for point sites or polyline sites,*L*_{−∞}(minimum),*L*_{0}(coordinate product),*L*_{∞}(maximum), any other*L*_{p}distance measures, the*Euclidean highway*and*Manhattan highway*distance measures (with faster travel on specified line segments), and the*triangular-grid*,*Karlsruhe*,*city*,*Köln*,*orbit-in*,*orbit-out*, and*spherical*distance measures for point sites,*angular-size*,*detour*, and*dilation*measures for polyline sites,*secant*,*catch distance*,*turn*,*left-turn*, and*Dubins path length*measures for half-line/ray/rooted-vector sites,- and weighted sums or products of these distance measures (including
*semi-Voronoi diagrams*),

- with or without distance contour lines (lines of equidistant points)
- in black and white or with colour-coded sites and regions.

Use Vorosketch wisely; it is not designed for speed or for accuracy. Vorosketch may miss features of the diagram that are narrower than a pixel, and I made no attempt to analyse or control the error margins in the distance computations. For standard Euclidean Voronoi diagrams of large numbers of points in the plane, AGB-DTVD and many other implementations are orders of magnitude faster and more precise.

Below you can find Vorosketch's source code, examples of input and output, some technical background on the different types of diagrams, and some answers to possibly-asked questions. As a teaser, below is a travel-time Voronoi diagram of 150 points, where travel along highways is ten times as fast as across the plane. For more large figures, see the gallery.

Vorosketch is written in C++, here is the source code. The current version is 0.24 beta of 17 September 2022. Vorosketch does not use any libraries apart from the standard template library, so it should be easy to compile with any C++ compiler. Please let me know if you have trouble compiling the code or if you find any bugs. Vorosketch is licensed under the Apache License, please see the source code for further details. To enable user-defined distances from userdistances.cpp (see the source code for an example), compile with the -D INCLUDE_USER_DISTANCES option.

Here are some examples of what Vorosketch can do. All examples were made with the “-r 1200” option to produce a drawing of 1200×1200 pixels in bmp-format. For convenience, we omit “-r 1200” from the command lines below. For publication on this webpage, the images were downsampled to 400×400 pixels and converted to png-format.

For a more complete list of options, please run Vorosketch with the -? option.

A standard Voronoi diagram of points in the Euclidean plane is produced by:

`vorosketch <sites >image.bmp`

Next: the same, with additive weights. The boundaries of the shaded areas are at equal weighted distance to their sites.

`vorosketch <sites >image.bmp`

Next: a power diagram (that is, using squared distances with additive weights). And why not add some colour, too, using Sasha Trubetskoy's beautiful colour scheme:

`vorosketch -m squared -c <sites >image.bmp`

With the -i and -g options one can add contour lines and shading:

`vorosketch -i 0.05 -c -g 0.2 <sites >image.bmp`

Euclidean distance with multiplicative weights:

`vorosketch -d <sites >image.bmp`

All you have seen so far, also works for polyline sites. Let us use colour (-c) and add the black region boundaries nonetheless (-b). As you can see, there are areas that are at equal distance from two sites:

`vorosketch -b -c <sites >image.bmp`

To cut the equal-distance areas in a sensible way, use the -s option to shorten the polylines slightly:

`vorosketch -b -c -s 0.01 <sites >image.bmp`

In the diagram below each site moves at constant speed in the direction indicated by the black “arm” that originates from its coloured dot. We define the distance from a point *q* to a site *p* as the distance one has to travel from *q* (at the same constant speed) to meet (catch) the moving site *p*. This leads to the following Voronoi diagram of moving points:

`vorosketch -m catch -i 0.025 -c -g 1.25 <sites >image.bmp`

We can also make a Voronoi diagram of points that move at different speeds, while the catcher moves at a fixed speed. For example:

Another interpretation of the same diagram is that a point *q* belongs to the site *p* that can reach *p* fastest, when the site's original speed and direction are subjected to an instantaneous directed acceleration of a fixed magnitude.

Now let the sites represent disks, oriented perpendicular to the drawing plane and perpendicular to the black arms drawn in the diagram. A disk at a point *p* emits light on the side of the black arm (whose end point we denote by *p'*) but not on the other side. Naturally, the light received by a point *q* in the plane is now inversely proportional with the secant of the angle *qpp'* and with the squared distance to *p*. We want to produce a diagram that shows, for each point in the plane, from which disk it receives most light. For this purpose, we use Vorosketch's ability to multiply distance measures:

`vorosketch -m squared*secant -i 0.05 -c -g 1 -b <sites >image.bmp`

Now suppose the sites are oriented line segments, and the distance to a point *p* is how much we would have to rotate the segment around its starting point in order to point towards *p*:

`vorosketch -m turn -b -c <sites >image.bmp`

For a vehicle, rotating can be an issue too. Here is a Voronoi diagram of 18 cars, each visualised as a rooted vector. The origin of the vector (marked by a coloured dot) marks the current location of the car, and the length indicates its minimum turning radius. The distance from a car to a point is defined as the length of the shortest route that the car can take to that point while only driving forward and subject to its minimum turning radius:

`vorosketch -m dubins -i 0.04 -b -c <sites >image.bmp`

In semi-Voronoi diagrams, distances are as usual, except that the distance between a point and a rooted-vector site is infinite if the point is not in front of it, that is, if the site would have to rotate more than 90 degrees to be directed at the point. Semi-Voronoi diagrams can be produced by multiplying the usual distance measure with the special “semi” distance, for example:
`vorosketch -m semi*euclidean -b -c <sites >image.bmp`

Some of the distance measures for point sites also work for segment or polyline sites. In addition, Vorosketch implements some distance measures that are defined specifically for segment and polyline sites. For exampe, suppose we say the closest object to *p* is the one that has the largest angular size as seen from *p*:

`vorosketch -m angle -b -c <sites >image.bmp`

Or the closest line segment *qr* to *p* is the one that minimises the detour if we would want to visit *p* on the way from *q* to *r*, that is, it minimises |*qp*| + |*pr*| - |*qr*|:

`vorosketch -m detour -i 0.1 -b -c <sites >image.bmp`

Under the Manhattan (L1) metric, the distance between two points is determined by the shortest path that consists only of vertical and horizontal segments. To make sure that regions that are equidistant to two sites are recognisable as such, use the -e option to catch rounding errors:

`vorosketch -m manhattan -e 0.01 -b -c <sites >image.bmp`

In fact, we can calculate Voronoi diagrams for any geometric *L*_{p} distance measure, even for *p* < 1, for example the minimum measure *L*_{−∞}, the product measure *L*_{0}, or *L*_{0.5}:

`vorosketch -m L-inf -b -c <sites >image.bmp`

`vorosketch -m L0 -b -c <sites >image.bmp`

`vorosketch -m L0.5 -b -c <sites >image.bmp`

Under the Karlsruhe metric, distances are determined by paths that consist only of radial segments (to/from the centre of the map) and circular arcs (with the circle centre in the centre of the map):

`vorosketch -m karlsruhe -e 0.01 -b -c -+ <sites >image.bmp`

In a real city, say, Köln, you would also have to take account that travelling in the city centre is slow. Suppose we can travel only on radial segments and circular arcs, with a speed that is linearly proportional to the distance from the centre, slightly faster on the circular arcs than on the radial axes (in effect, this turns distances into L1-distances in polar coordinate space with a logarithmic radius axis, with the actual centre point becoming unreachable). We use Vorosketch's ability to add up multiple distance measures:

`vorosketch -m 0.9*azimuth+logradius -b -c -e 0.01 <sites >image.bmp`

Now suppose the closest site to *p* is the one that would require least energy to launch a space ship that reaches *p* from that site, subject to gravitation towards the centre:

`vorosketch -m orbitout -b -c -+ <sites >image.bmp`

The closest site to *p* is the one that would require least energy to launch a space ship that reaches the site from *p*:

`vorosketch -m orbitin -b -c -+ <sites >image.bmp`

Very different yet is the diagram that results, if one can travel in all directions, but travel on highways is faster than cross-country. For this diagram, the -i option is used to show contour lines at intervals of 0.05:

`vorosketch -m highway -b -c -e 0.01 -i 0.05 -g 0.5 <sites >image.bmp`

Finally, a second-order Voronoi diagram:

`vorosketch -2 -b -c <sites >image.bmp`

Vorosketch simply computes, for each pixel, the distance to each site. Of course that is not a particularly efficient way to compute the diagrams, but it works for any type of diagram, and for small diagrams to be used as illustrations, it is fast enough. If there are multiple closest sites, Vorosketch records *two* of them (if there are more than two, Vorosketches takes the first two by order in the input file). The reason that only the two closest sites are considered, is that so far, that was enough for my purposes. To handle situations in which there are more than two closest sites, we would first have to decide how we would want that to be visualised. I am open to suggestions if there is a need for it.

The case of Euclidean distance is well-known. Contour lines of point sites are concentric circles. Bisectors of point sites are straight lines and/or hyperboles (in the additively/subtractively weighted case) or Apollonian circles (in the multiplicatively/divisively weighted case). With unweighted line segment sites, contour lines consist of concentric circular arcs and straight line segments that advance at unit speed. Bisectors that are traced out by advancing circular arcs meeting advancing straight line segments take the form of parabolic curves.

The *catch* distance from a point *q* to a ray from *p* through *p'* is defined as the secant of the angle *qpp'* times half of the distance |*qp*|, provided the angle *qpp'* is positive (otherwise the catch distance is infinite). Contours are circles with *p* on the boundary and the centre points on the ray. The distance models two very different “semantics”. First (as is obvious from the contours): how long does it take for somebody that starts at *q* and can move in any direction to catch up with somebody that starts at *p* and moves on the ray through *p'*, if both move at the same speed? Second (this is how I came up with the definition of this distance measure): if the ray from *p* models a small plate of some small standard width ε, with normal oriented towards *p'*, then what is the reciprocal of the angular size of the front (the side oriented towards *p'*) of the plate in the view from *q*? In other words, how (in)visible is the plate from *q*?
When |*pp'*| = 1, we can rewrite the catch distance as ( (*q*−*p*).(*q*−*p*) ) / ( (*q*−*p*).(*p'*−*p*) ). Two sites *a* and *b* are at equal distance from *q* when ( (*q*−*a*).(*q*−*a*) ) * ( (*q*−*b*).(*b'*−*b*) ) − ( (*q*−*b*).(*q*−*b*) ) * ( (*q*−*a*).(*a'*−*a*) ) = 0 and the secants are positive. Thus, the bisector of *a* and *b* is (part of) the zero set of a third-degree polynomial.

With the *turn* distance, the distance between a point *P* and a ray from *A* through *B* is the angle between *AP* and *AB*. If the rays are not weighted, then the bisectors consist of straight line segments, circular arcs (see Alegría et al.) and hyperbolic segments (see Haverkort and Klein).

With the *Dubins* distance, a site is a vector *v* rooted at a point *p*, representing a car at *p*, oriented in the direction *v*, with minimum turning radius |*v*|. The distance from such a site to a point *q* in the plane is the length of the shortest path that the car can drive from *p* to *q* without driving backwards. Depending on the location of *q* relative to *p* and *v*, the path consists of up to two circular arcs and one straight segment—see Bui and Boissonnat's 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.

With the angle (or: angular-size) distance measure, a contour line is the locus of points from which a site looks equally large in the field of view. If the sites are line segments, then, by the inscribed-angle theorem, the contour lines are (non-concentric) circles with the segment endpoints on the boundary. The bisectors are interesting curves with inflection points, whose properties I have not investigated yet.

With the detour and dilation distance measures, the “distance” between a point *X* and a line segment *AB* is defined as |*AX*|+|*XB*|-|*AB*|; in the case of dilation, this number is subsequently divided by |*AB*|. Thus, a contour line (locus of points equidistant to *AB*) consists of the points *X* with constant sum |*AX*|+|*XB*|, which constitutes an ellipse. I have not thought about what shapes the bisectors have yet; any thoughts welcome!.

Klein and I found that it makes sense to define the *geometric L*_{0} distance between two points as the area of their axis-alligned bounding box (or in other words: the product of their distances on the coordinate axes). If the sites are not weighted, we found that bisectors consist of hyperbolic curves and straight lines.

The Voronoi diagram under the Karlsruhe metric (also known as Moscow or Amsterdam metric) has been discussed by Klein. The city, Köln, orbit-in and orbit-out metrics are, to the best of my knowledge, novelties that I added for fun.

The city and Köln metrics are simply the Euclidean and the Manhattan metrics after the coordinates have been converted to a polar coordinate system with a logarithmic radius axis, and where the azimuth axis wraps around from 2π to 0. If the sites are not weighted, the resulting bisectors consist of straight line segments in the logarithmic polar coordinate system. Thus, in the original coordinate system (that is, what you see), the bisectors consist of segments of straight lines through the origin and of circles and logarithmic spirals around the origin.

The orbit distance measures are based on the laws of physics. (Of course this does not mean the setting is realistic in any way—in particular, in space, the sites themselves would also be orbiting the centre, they would not be stationary.)
The speed of a small object at position *X* in an elliptic orbit around a large object in the origin *O* is proportional to √(1/|*OX*| − 1/(2*a*)), where *a* is the length of the semi-major axis of the ellipse, and therefore the kinetic energy is proportional to 1/|*OX*| − 1/(2*a*). If we want to launch a space ship from *X* such that it will reach *Y*, we have to give it a direction and a kinetic energy that puts it into an orbit that reaches *Y*. Since this energy increases with *a*, we are looking for the orbit that includes *X* and *Y*, has *O* as one of its focal points, and has the smallest possible semi-major axis (non-elliptical, that is, open-ended orbits do not need to be considered: they would require more energy in any case). Note that, in an ellipse, *a* is exactly 1/4 of the length of the path *OXFYO*, where *F* is the second focal point. Therefore, to minimise *a*, we minimise the length of *OXFYO*, which we achieve by putting *F* on the line segment *XY*. Thus, the required energy is proportional to 1/|*OX*| − 2/(|*OX*|+|*XY*|+|*YO*|), which is straightforward to compute.

With the orbit-out measure, *X* is an input site, for which the first term is fixed, and a contour line (locus of points at equal distance to this site) consists of points *Y* with constant |*XY*|+|*YO*|: that is an ellipse with focal points *X* and *O*.
With the orbit-in measure, *Y* is the input site, and the situation seems to be more complicated.
In either case, I have not given the shape of the resulting bisectors any thought yet, and I have not figured out if the distance measures satisfy the triangle inequality yet; I would appreciate any thoughts.

With unweighted point sites in the highway distance, finding the shortest path between points becomes part of the challenge. Each segment of such a shortest path has two endpoints. Some of these endpoints are inherent to the input configuration (point sites, endpoints of highway segments, and intersection points of highway segments), the others (where one enters or leaves a highway at an arbitrary point) we may call *free*. Two observations help us limit the complexity of finding shortest paths.

First, when we enter or leave a highway with speed *s*>1 at a free point, we always enter or leave it under an angle with cosine equal to 1 divided by *s*. (Here, the speed of travelling cross-country is taken to be 1.) One can easily verify that this angle is optimal: if we would enter the highway at a bigger angle, we would be making too much of a detour; if we would enter the highway at a smaller angle, we would not be taking advantage of the highway enough.

Second, to find shortest paths, we do not need to consider cross-country segments with two free endpoints. To see this, consider any cross-country segment *e* with two free endpoints *d* and *f*, that is, a segment that shortcuts from one highway segment *C* to another highway segment *G*. Let *c*=*bd* on *C* be the previous segment of the shortest path, and let *g*=*fh* be the next segment of the shortest path. Now suppose we move *e* to the left or to the right while maintaining its slope and maintaining the connection to *C* (by moving *d* on *C*) and to *G* (by moving *f* on *G*), so that we maintain the angles on which we leave *C* and enter *G* as per the first observation. Then the total path length changes linearly. Thus, in at least one direction, the total path length will not increase and we can move *e* until one of its endpoints, say *d* (the case of *f* is symmetric), hits either
(i) a site, a highway endpoint or an intersection point *p*, or
(ii) the free endpoint *b* that started *c*.
In case (i), the other endpoint of *e*, that is, *f* on *G*, must be the free endpoint where a ray from *p* hits *G* at the correct angle. Given *p* and *G*, there are at most two such points, so we can actually precompute all such candidate free endpoints.
In case (ii), *c* is now eliminated, and the previous segment of the shortest path, which is also a cross-country segment, must be collinear with *e*. Thus, we can treat the previous segment and *e* together as a single shortest path segment and continue our adaptations of the shortest path with two segments less than before. In the end, we can transform any shortest path into a shortest path in which each cross-country segment has at least one non-free endpoint, and all free endpoints come from a finite set that can be computed easily.

These observations allow us to precompute much of the shortest paths that we need to generate the diagram (see El Shawi et al. for further details). After that, all that is left to do is figuring out, for each pixel, what is the first input point or highway segment on that path. To do this efficiently, Vorosketch uses a heuristic solution: it first computes the eligible input points and highway segments for each cell in a 100×100 grid (using conservative bounds on the distances that can be computed fast); then, for each pixel in the cell, only shortest paths that go through these input points or highway segments need to be considered.

Contours of a point site are circular arcs and straight line segments advancing at unit speed, where the straight line segments originate from highway sections at an angle with sine equal to 1 divided by the highway speed. However, the circular arcs in a contour may have different radius, as they may grow from given points (for example, highway endpoints) that are at different distances to the starting point (site). Thus, bisectors can be straight line segments (as in unweighted Voronoi diagrams), hyperbolic arcs (as in additively weighted Voronoi diagrams) or parabolic arcs (as in line segment Voronoi diagrams).

For more on highway Voronoi diagrams, one may consult Bae and Chwa.

**What is new since version 0.20?**

- Since version 0.24, Vorosketch supports distances on the sphere rendered according to an equal-area azimuthal, cylindrical, or elliptical projection.
- Since version 0.24, distance in a triangular grid is supported.
- Since version 0.24, colour schemes can be read from a file.
- Since version 0.24, distances can be shown by a customisable colour gradient (heat map).
- Since version 0.23, the
*Manhattan highway*distance is supported, in which non-highway travel is measured by the L1 metric. - Since version 0.22, for most distance measures Vorosketch makes a pre-scan of a grid of squares, in which it obtains bounds on the distances to the sites from each square. These bounds are then used to determine for which pixels and sites it makes sense to compute the distance exactly, and for which it does not. Depending on the distance measure and the input, this may speed up Vorosketch by an order of magnitude or more.
- Since version 0.22, distance functions can be composed on the command line from built-in distance functions and a number of standard mathematical operators. This should also work with the highway distance now.
- Since version 0.22, there are many new colour options, including a choice of palettes, heuristics (experimental: subject to change) to assign colours as to maximise or minimise contrast between adjacent regions, shading based on the gradient of the distance function, and saturation based on the number of neighbours of a region.
- Since version 0.22, Vorosketch can generate random sites (experimental: subject to change).
- Sinve version 0.22, the unit speed
*v*for the catch distance with mixed-speed sites is no longer selected by “-m mixedcatch -i*v*”, but by “-m catch@*v*”.

**The output seems wrong. Why?**

Sometimes I find bugs and fix them, or I adapt the user interface to catch problematic combinations of parameter settings. If you discover a bug, please check if you are using the newest version of Vorosketch. If you compiled with the -ffast-math option, then try without: Vorosketch relies on correct handling of ∞, −∞ and not-a-number. If the bug is still there, please let me know.

**Why does shading by distance gradient (the -z option) does not work?**

It should work. However, the distances are not scaled automatically. Depending on the input, some distance measures tend to produce rather small values for points in the [−1,1]x[−1,1] square, so that the distance gradients are small and the distance landscape appears flat. In particular, the squared Euclidean, secant, and catch distances may have this problem (whereas the turn and leftturn distances tend to have too much variation, producing steep gradients). To remedy this, scale the distances. For example, use “-m squared*5” instead of “-m squared”.

**Could diagram X not be computed faster?**

Certainly. The computation of the angular-size, detour, dilation, catch (with sites moving at different speeds), and Dubins Voronoi Diagrams is particularly slow, because good bounds for the pre-scan heuristic have not been implemented. Since I have to budget the time I spent programming, I try to be selective in what computations I choose to optimise. So far, most diagrams that I wanted to produce were produced within seconds. Theoretically, speeding it up further only makes sense for types of diagrams of which hundreds of thousands are going to be produced. If you are going to produce thousands of Voronoi diagrams that need to be rendered faster, please let me know and I might be tempted to speed up my efforts to speed it up.

**Why do I get an error: “Distance function is not compatible with or has not been implemented for any type of site”?**

This happens if one combines distances whose domains (types of sites for which they are defined and implemented) are incompatible. For example, using “-m karlsruhe+dilation” results in this error, as the Karlsruhe distance has only been implemented for point sites, whereas the dilation distance is only defined for polyline sites.

**Why do I get an error: “missing distance name”?**

This happens if, in the argument of the -m option, a distance name or formula is missing where one is expected. For example, using “-m euclidean+” results in this error, because the second operand of the + operator is missing.

**Why do I get an error: “negative weight detected in input”?**

If you use divisive (multiplicative) weighting (the -d option), weights are not allowed to be negative. The reason is that lower bounds on distances would become upper bounds and vice versa: the implementation would need to handle this and currently it does not. Please let me know if you had a use for negative divisive (multiplicative) weights.

**Why do I get an error: “negative value -1.0842e-19 for upper bound on first operand to power operator detected” (or a similar message)?**

The power operator does not accept negative values for the first operand (for taking squares, use the sq function instead). If you suspect that the negative value is the result of a rounding error, use the maximum operator | to round negative values up to 0, for example: -m “(ln(karlsruhe/euclidean)|0)^1.5”

**Why do I get an error: “too many sites”?**

Vorosketch is designed for relatively small inputs, it is not designed to scale well. The code sets a somewhat arbitrary maximum of 20,000 sites (less when using highway distance). Not very far above that, Vorosketch would crash because it uses a quadratic amount of memory (which is entirely avoidable, but would require some programming effort). With some inputs, distance measures and options, Vorosketch would be impossibly slow with far fewer sites already, due to the arithmetic operations involved.

**Why is the output in this inefficient bmp-format?**

Because I had C++ code for that lying around already and I did not get round to studying PNG yet. This is somewhere on my to-do list, but all the way at the bottom, to be honest.

**Mr. Haverkort, could you please add feature X?**

Just ask me and I will have a look.

vorosketch.txt · Last modified: 2022/10/10 21:00 by administrator