Category Archives: Julia

SyntheticGrids.jl: Part 1

By: Eric Perim

Re-posted from: https://invenia.github.io/blog/2018/06/27/syntheticgrids-part-1/

Background

For the package repository, visit Github.

It should come as no surprise that electricity plays a vital role in many aspects of modern life. From reading this article, to running essential hospital equipment, or powering your brand-new Tesla, many things that we take for granted would not be possible without the generation and transmission of electrical power. This is only possible due to extensive power grids, which connect power producers with consumers through a very complex network of towers, transmission lines, transformers etc. Needless to say, it is important to understand the peculiarities of these systems in order to avoid large scale blackouts, or your toaster burning out due to a fluctuation in the current.

Power grid research requires testing in realistic, large-scale, electric networks. Real power grids may have tens of thousands of nodes (also called buses), interconnected by multiple power lines each, spanning hundreds of thousands of square kilometers. In light of security concerns, most information on these power grids is considered sensitive and is not available to the general public or to most researchers. This has led to most power transmission studies being done using only a few publicly available test grids 1, 2. These test grids tend to be too small to capture the complexity of real grids, severely limiting the practical applications of such research. With this in mind, there has recently been an effort in developing methods for building realistic synthetic grids, based only on publicly available information. These synthetic grids are based on real power grids and present analogous statistical properties—such as the geographic distribution of load and generation, total load, and generator types—while not actually exposing potentially sensitive information about a real grid.

The pioneers in treating power grids as networks were Watts and Strogatz3, when they pointed out that electric grids share similarities with small-world networks: networks that are highly clustered, but exhibit small characteristic path lengths due to a few individual nodes being directly connected to distant nodes (see Figure 1). This type of network is very useful in explaining social networks— see six degrees of Kevin Bacon—but, despite similarities, power grids differ from small-world networks 4,5. If you are looking for an extensive list of studies on power grids, Pagani and Aiello 6 is a good place to start.

Network topologies
Some examples of different network topologies containing 20 nodes and 40 edges. (a) Small-world; (b) random; (c) scale-free (exponent 2). For more details, see Watts and Strogatz3.

In order to study the dynamic properties of electric grids, some research has adopted simplified topologies, such as tree structures 7 or ring structures 8, which may fail to capture relevant aspects of the system. Efforts to build complete and realistic synthetic grids are a much more recent phenomenon. The effort of two teams is particularly relevant for this post, namely, Overbye’s team 9,10,11 and Soltan and Zussman 12.

Considering the potential impact of synthetic grids in the study of power grids and the recency of these approaches, we at Invenia Labs have developed SyntheticGrids.jl, an open source Julia package. The central idea of SyntheticGrids.jl is to provide a standalone and easily expandable framework for generating synthetic grids, adopting assumptions based on research by Overbye’s and Zussman’s teams. Currently, it only works for grids within the territory of the US, but it should be easily extendable to other regions, provided there is similar data available.

There are two key sources of data for the placement of loads and generators: USA census data and EIA generator survey data. The former is used to locate and size loads, while the latter is used for generators. Since there is no sufficiently granular location-based consumption data available, loads are built based on population patterns. Load has a nearly linear correlation with population size 12, so we adopt census population as a proxy for load. Further, loads are sited at each zip code location available in the census data. When placing generators, the EIA data provides us with all the necessary information, including geographic location, nameplate capacity, technology type, etc. This procedure is completely deterministic, since we want be as true as possible to the real grid structure, i.e. we want to use an unaltered version of the real data.

Our package treats power grids as a collection of buses connected by transmission lines. Buses can be either load or generation buses. Each generation bus represents a different power plant, so it may group several distinct generators together. Further, buses can be combined into substations, providing a coarse-grained description of the grid.

The coarse-graining of the buses into substations, if desired, is done via a simple hierarchical clustering procedure, as proposed by Birchfield et al. 9. This stochastic approach starts with each bus being its own cluster. At each step, the two most similar clusters (determined by the similarity measure of choice) are fused into one, and these steps continue until a stopping criterion has been reached. This allows the grouping of multiple load and generator units, similarly to what is actually done by Independent System Operators (ISOs).

In contrast to loads and generators, there is no publicly available data on transmission lines, so we have to adopt heuristics. The procedure implemented in the package is based on that proposed by Soltan and Zussman 12. It adopts several realistic considerations in order to stochastically generate the whole transmission network, which are summarised in the following three main principles:

The degree distributions of power grids are very similar to those of scale-free networks [see: Scale-free network], but grids have less degree 1 and 2 nodes and do not have very high degree nodes.

It is inefficient and unsafe for the power grids to include very long lines.

Nodes in denser areas are more likely to have higher degree.

Currently, SyntheticGrids.jl allows its generated grids to be directly exported to pandapower, a Python-based powerflow package. Soon, an interface with PowerModels.jl, a Julia-based powerflow package, will also be provided.

In the second part we will go over how to use the main features of the package.

References

  1. Power Systems Test Case Archive (UWEE) – http://www2.ee.washington.edu/research/pstca/ 

  2. Power Cases – Illinois Center for a Smarter Electric Grid (ICSEG) – http://icseg.iti.illinois.edu/power-cases 

  3. Watts, D. J., & Strogatz, S. H. (1998). Collective dynamics of ‘small-world’ networks. nature, 393(6684), 440-442. Chicago  2

  4. Hines, P., Blumsack, S., Sanchez, E. C., & Barrows, C. (2010, January). The topological and electrical structure of power grids. In System Sciences (HICSS), 2010 43rd Hawaii International Conference on (pp. 1-10). IEEE. 

  5. Cotilla-Sanchez, E., Hines, P. D., Barrows, C., & Blumsack, S. (2012). Comparing the topological and electrical structure of the North American electric power infrastructure. IEEE Systems Journal, 6(4), 616-626. 

  6. Pagani, G. A., & Aiello, M. (2013). The power grid as a complex network: a survey. Physica A: Statistical Mechanics and its Applications, 392(11), 2688-2700. 

  7. Carreras, B. A., Lynch, V. E., Dobson, I., & Newman, D. E. (2002). Critical points and transitions in an electric power transmission model for cascading failure blackouts. Chaos: An interdisciplinary journal of nonlinear science, 12(4), 985-994. 

  8. Parashar, M., Thorp, J. S., & Seyler, C. E. (2004). Continuum modeling of electromechanical dynamics in large-scale power systems. IEEE Transactions on Circuits and Systems I: Regular Papers, 51(9), 1848-1858. 

  9. Birchfield, A. B., Xu, T., Gegner, K. M., Shetye, K. S., & Overbye, T. J. (2017). Grid structural characteristics as validation criteria for synthetic networks. IEEE Transactions on power systems, 32(4), 3258-3265. Chicago  2

  10. Birchfield, A. B., Gegner, K. M., Xu, T., Shetye, K. S., & Overbye, T. J. (2017). Statistical considerations in the creation of realistic synthetic power grids for geomagnetic disturbance studies. IEEE Transactions on Power Systems, 32(2), 1502-1510. Chicago 

  11. Gegner, K. M., Birchfield, A. B., Xu, T., Shetye, K. S., & Overbye, T. J. (2016, February). A methodology for the creation of geographically realistic synthetic power flow models. In Power and Energy Conference at Illinois (PECI), 2016 IEEE (pp. 1-6). IEEE. 

  12. Soltan, Saleh, and Gil Zussman. “Generation of synthetic spatially embedded power grid networks.” arXiv:1508.04447 [cs.SY], Aug. 2015.  2 3

SyntheticGrids.jl: Part 2

By: Eric Perim

Re-posted from: https://invenia.github.io/blog/2018/06/27/syntheticgrids-part-2/

Usage

For the package repository, visit Github.

In the first part, we discussed the motivation and model behind SyntheticGrids.jl. In this post we show how to use it.

To use SyntheticGrids.jl, Julia 0.6.1 or newer is required. Once Julia is properly installed, the package can be installed via

julia> Pkg.add("SyntheticGrids")

This should take care of all dependencies. In order to check if the package has been properly installed, use

julia> Pkg.test("SyntheticGrids")

A (very) simple test example

As an introduction to the package, we start by automatically generating a small, but complete grid.

julia> using SyntheticGrids

julia> grid = Grid(false);

This command generates a complete grid corresponding to the region contained in the box defined by latitude [33, 35] and longitude [-95, -93] (default values). It automatically places loads and generators and builds the transmission line network (we will soon see how to do each of these steps manually). Here, false determines that substations will not be created. Note the addition of the semicolon, ;, at the end of the command. This has just cosmetic effect in suppressing the printing of the resulting object in the REPL. Even a small grid object corresponds to a reasonably large amount of data.

A Grid object has several attributes that can be inspected. First, let’s look at the buses:

julia> length(buses(grid))
137

julia> buses(grid)[1]
LoadBus(
	id=1,
	coords=LatLon(lat=33.71503°, lon=-93.166445°),
	load=0.17400000000000002
	voltage=200,
	population=87,
	connected_to=Set{Bus}(...)
	connections=Set{TransLine}(...)
)

julia> buses(grid)[end]
GenBus(
	id=137
	coords=LatLon(lat=34.4425°, lon=-93.0262°),
	generation=56.0
	voltage=Real[115.0],
	tech_type=AbstractString["Conventional Hydroelectric"],
	connected_to=Set{Bus}(...)
	connections=Set{TransLine}(...)
	pfactor=0.9
	summgen=61.8
	wintgen=62.0
	gens=SyntheticGrids.Generator[SyntheticGrids.Generator(LatLon(lat=34.4425°, lon=-93.0262°), Real[115.0], "Conventional Hydroelectric", 28.0, 0.9, 15.0, 30.9, 31.0, "1H", "OP"), SyntheticGrids.Generator(LatLon(lat=34.4425°, lon=-93.0262°), Real[115.0], "Conventional Hydroelectric", 28.0, 0.9, 15.0, 30.9, 31.0, "1H", "OP")]
)

We see that our grid has a total of 137 buses (see Figure 2 for a visualisation of the result). The first is a load bus (LoadBus). The values of the attributes connected_to and connections are not explicitly printed. However, the printing of (...) indicates that those sets have been populated (otherwise, they would be printed as ()).

Synthetic grids
Visualisation of two grids generated using the procedure described here. Notice that both present the same bus locations, as their placement is entirely deterministic. The transmission line topology however is different in each case, as it is generated through an stochastic process. Note that the generated grids are non-planar.

The last bus of the list corresponds to a generator (GenBus). One important thing to notice here is that it contains an attribute called gens, which is an array of Generator-type objects. GenBuses represent power plants, which may (or may not, as is the case here) contain several different generating units. These individual generating units are stored within the gens attribute.

We can also inspect the transmission lines:

julia> length(trans_lines(grid))
167

julia> trans_lines(grid)[1]
TransLine(
	connecting: (LoadBus(
	id=3,
	coords=LatLon(lat=33.889332°, lon=-93.097793°),
	load=8.18
	voltage=100,
	population=4090,
	connected_to=Set{Bus}(...)
	connections=Set{TransLine}(...)
), LoadBus(
	id=1,
	coords=LatLon(lat=33.71503°, lon=-93.166445°),
	load=0.17400000000000002
	voltage=200,
	population=87,
	connected_to=Set{Bus}(...)
	connections=Set{TransLine}(...)
)),
	impedance=0.9175166312451004,
	capacity=1400
)

There are 167 transmission lines in our grid. By looking at the first one, we see that they are defined by a tuple of Bus-type objects (here both are LoadBuses), by an impedance value (here taken as Real, since the package has been developed with DC OPF in mind), and a current carrying capacity value.

The adjacency matrix of the system can also be easily accessed:

julia> adjacency(grid)
137×137 SparseMatrixCSC{Bool,Int64} with 334 stored entries:
  [3  ,   1]  =  true
  [6  ,   1]  =  true
  [15 ,   1]  =  true
  [34 ,   1]  =  true
  [35 ,   1]  =  true
  [4  ,   2]  =  true
  
  [54 , 135]  =  true
  [58 , 135]  =  true
  [67 , 135]  =  true
  [73 , 136]  =  true
  [42 , 137]  =  true
  [46 , 137]  =  true

Notice that we use a sparse matrix representation for better efficiency.

Substations can also be inspected, but we did not create any, so the result should be empty:

julia> substations(grid)
0-element Array{SyntheticGrids.Substation,1}

That can be remedied by changing the boolean value when creating the grid:

julia> grid = Grid(true);

julia> length(substations(grid))
43

julia> substations(grid)[end]
Substation(
	id=43
	coords=LatLon(lat=34.412130070351765°, lon=-93.11856562311557°),
	voltages=Real[115.0],
	load=0,
	generation=199.0,
	population=0,
	connected_to=Set{Substation}(...)
	grouping=SyntheticGrids.Bus[GenBus(
	id=137
	coords=LatLon(lat=34.4425°, lon=-93.0262°),
	generation=56.0
	voltage=Real[115.0],
	tech_type=AbstractString["Conventional Hydroelectric"],
	connected_to=Set{Bus}(...)
	connections=Set{TransLine}(...)
	pfactor=0.9
	summgen=61.8
	wintgen=62.0
	gens=SyntheticGrids.Generator[SyntheticGrids.Generator(LatLon(lat=34.4425°, lon=-93.0262°), Real[115.0], "Conventional Hydroelectric", 28.0, 0.9, 15.0, 30.9, 31.0, "1H", "OP"), SyntheticGrids.Generator(LatLon(lat=34.4425°, lon=-93.0262°), Real[115.0], "Conventional Hydroelectric", 28.0, 0.9, 15.0, 30.9, 31.0, "1H", "OP")]
), GenBus(
	id=135
	coords=LatLon(lat=34.570984°, lon=-93.194425°),
	generation=75.0
	voltage=Real[115.0],
	tech_type=AbstractString["Conventional Hydroelectric"],
	connected_to=Set{Bus}(...)
	connections=Set{TransLine}(...)
	pfactor=0.9
	summgen=75.0
	wintgen=75.0
	gens=SyntheticGrids.Generator[SyntheticGrids.Generator(LatLon(lat=34.570984°, lon=-93.194425°), Real[115.0], "Conventional Hydroelectric", 37.5, 0.9, 20.0, 37.5, 37.5, "10M", "OP"), SyntheticGrids.Generator(LatLon(lat=34.570984°, lon=-93.194425°), Real[115.0], "Conventional Hydroelectric", 37.5, 0.9, 20.0, 37.5, 37.5, "10M", "OP")]
), GenBus(
	id=136
	coords=LatLon(lat=34.211913°, lon=-93.110963°),
	generation=68.0
	voltage=Real[115.0],
	tech_type=AbstractString["Hydroelectric Pumped Storage", "Conventional Hydroelectric"],
	connected_to=Set{Bus}(...)
	connections=Set{TransLine}(...)
	pfactor=0.95
	summgen=68.0
	wintgen=68.0
	gens=SyntheticGrids.Generator[SyntheticGrids.Generator(LatLon(lat=34.211913°, lon=-93.110963°), Real[115.0], "Conventional Hydroelectric", 40.0, 0.95, 15.0, 40.0, 40.0, "1H", "OP"), SyntheticGrids.Generator(LatLon(lat=34.211913°, lon=-93.110963°), Real[115.0], "Hydroelectric Pumped Storage", 28.0, 0.95, 15.0, 28.0, 28.0, "1H", "OP")]
)]
)

By changing the boolean value to true we now create substations (with default values; more into that later) and can inspect them.

A more complete workflow

Let’s now build a grid step by step. First, we start by generating an empty grid:

julia> using SyntheticGrids

julia> grid = Grid()
SyntheticGrids.Grid(2872812514497267479, SyntheticGrids.Bus[], SyntheticGrids.TransLine[], SyntheticGrids.Substation[], Array{Bool}(0,0), Array{Int64}(0,0))

Notice that one of the attributes has been automatically initialised. That corresponds to the seed which will be used for all stochastic steps. Control over the seed value gives us control over reproducibility. Conversely, that value could have been specified via grid = Grid(seed).

Now let’s place the load buses. We could do this by specifying latitude and longitude limits (e.g.: place_loads_from_zips!(grid; latlim = (30, 35), longlim = (-99, -90))), but let’s look at a more general way of doing this. We can define any function that receives a tuple containing a latitude–longitude pair and returns true if within the desired region and false otherwise:

julia> my_region(x::Tuple{Float64, Float64}, r::Float64) = ((x[1] - 33)^2 + (x[2] + 95)^2 < r^2)
my_region (generic function with 1 method)

julia> f(x) = my_region(x, 5.)
f (generic function with 1 method)

julia> place_loads_from_zips!(grid, f)

julia> length(buses(grid))
3287

Here, my_region defines a circle (in latitude-longitude space) of radius r around the point (33, -95). Any zip code within that region is added to the grid (to a total of 3287) as a load bus. The same can be done for the generators:

julia> place_gens_from_data!(grid, f)

julia> length(buses(grid))
3729

This command adds all generators within the same region, bringing the total amount of buses to 3729.

We can also manually add extra load or generation buses if we wish:

julia> a_bus = LoadBus((22., -95.), 12., 200, 12345)
LoadBus(
	id=-1,
	coords=LatLon(lat=22.0°, lon=-95.0°),
	load=12.0
	voltage=200,
	population=12345,
	connected_to=Set{Bus}()
	connections=Set{Transline}()
)

julia> SyntheticGrids.add_bus!(grid, a_bus)

julia> length(buses(grid))
3730

The same works for GenBuses.

Once all buses are in place, it is time to connect them with transmission lines. This can be done via a single function (this step can take some time for larger grids):

julia> connect!(grid)

julia> length(trans_lines(grid))
0

This function goes through the stochastic process of creating the system’s adjacency matrix, but it does not create the actual TransLine objects (hence the zero length). That is done via the create_lines! function. Also note that connect! has several parameters for which we adopted default values. For a description of those, see ? connect.

Before we create the lines, it is interesting to revisit adding new buses. Now that we have created the adjacency matrix for the network, we have two options when adding a new bus: either we redo the connect! step in order to incorporate the new bus in the grid, or we simply extend the adjacency matrix to include the new bus (which won’t have any connections). This is controlled by the reconnect keyword argument that can be passed to add_bus!. In the former case, one uses reconnect = false (the default option); connections can always be manually added by editing the adjacency matrix (and the connected_to fields of the involved buses).

Once the adjacency matrix is ready, TransLine objects are created by invoking the create_lines! function:

julia> SyntheticGrids.create_lines!(grid)

julia> length(trans_lines(grid))
4551

We have generated the connection topology with transmission line objects. Finally, we may want to coarse-grain the grid. This is done via the cluster! function, which receives as arguments the number of each type of cluster: load, both load and generation or pure generation. This step may also take a little while for large grids.

julia> length(substations(grid))
0

julia> cluster!(grid, 1500, 20, 200)

julia> length(substations(grid))
1700

At this point, the whole grid has been generated. If you wish to save it, the functions save and load_grid are available. Please note that the floating-point representation of numbers may lead to infinitesimal changes to the values when saving and reloading a grid. Besides precision issues, they should be equivalent.

julia> save(grid, "./test_grid.json")

julia> grid2 = load_grid("./test_grid.json")

Some simple statistics can be computed over the grid, such as the average node degree and the clustering coefficient:

julia> mean_node_deg(adjacency(grid))
2.4402144772117964

julia> cluster_coeff(adjacency(grid))
0.08598360707539486

The generated grid can easily be exported to pandapower in order to carry out powerflow studies. The option to export to PowerModels.jl should be added soon.

julia> pgrid = to_pandapower(grid)
PyObject This pandapower network includes the following parameter tables:
   - load (3288 elements)
   - trafo (913 elements)
   - ext_grid (1 elements)
   - bus_geodata (3730 elements)
   - bus (3730 elements)
   - line (3638 elements)
   - gen (1397 elements)

Conclusion

Hopefully, this post helped as a first introduction to the SyntheticGrids.jl package. There are more functions which have not been mentioned here; the interested reader should refer to the full documentation for a complete list of methods. This is an ongoing project, and, as such, several changes and additions might still happen. The most up-to-date version can always be found at Github.

Bezier Moi

By: cormullion

Re-posted from: https://cormullion.github.io/blog/2018/06/21/bezier.html

This post was inspired mostly by the arrival of the new Julia package Literate.jl by Fredrik Ekre. Literate lets you write a Julia source file, a Markdown blog post, and a Jupyter notebook all at the same time. It’s magic, or, at least, indistinguishable from magic… As a result, you might be able to find a Jupyter notebook version of this Markdown-converted Julia source file somewhere nearby (look in the github repo for this blog, perhaps). The code uses the following packages, Literate, Luxor, Colors, Roots, Fontconfig, DataFrames, Iterators, ColorSchemes, and should work in Julia v0.6 (some packages such as Colors have yet to be updated to work with version 0.7).

Bézier moi!

Luxor provides some support for Bézier curves, but there’s no room for any more documentation—it’s already too big. So the intention of this post is to provide some of the missing information about what’s up with Bézier curves and how you might use them. And I confess in advance to wasting some of your bandwidth with some pointlessly colourful graphics.

There’s not a lot of mathematical material here, but fortunately the internet is awash with high-quality information about Bézier curves. The two articles you should definitely read instead of this post, or at least before, are:

Historical background

The story is quite well known: two engineers employed in the French car industry, Paul de Faget de Casteljau, at Citroen, and Pierre Etienne Bézier, at Renault, worked—mostly independently—on the mathematics of curves in the early 1960s, as the industry made its first tentative steps towards using computers for design and production.

Citroen DS

(This is the Citroen DS, which was once voted the most beautiful car ever made, apparently. Image from the WikiMedia Commons photographed by Klugschnacker)

De Casteljau and Bézier were interested in mathematical tools that would allow designers to intuitively construct and manipulate complex shapes. This problem was especially critical for “free–form” shapes that couldn’t easily be specified by centre points, axes, angles, and dimensions. The motivation was also partly to replace the laborious, variable, and expensive process of sculpting clay models to specify the desired shape.

De Casteljau found some resistance when his mathematical researches were introduced into the design studio. He observed that:

the designers were astonished and scandalized. Was it some kind of joke? It was considered nonsense to represent a car body mathematically. It was enough to please the eye, the word ‘accuracy’ had no meaning. [quoted in Farouki]

Bézier popularized, but did not actually create, what we know today as the Bézier curve. He mainly developed the notation, and devised the idea of nodes with attached “control handles”, which the designers could use to adjust the shapes as easily as they used the turn indicators on their beloved Citroen DSs.

De Casteljau is also remembered for the algorithm that bears his name.

The shortest distance?

With the lingering thought of old Renaults and Citroens in mind, it’s tempting to think of a Bézier curve as a line that takes, not the shortest distance between two points, but takes instead the scenic route.

Let’s define four points:

using Luxor
P1  = Point(-200, 0)
CP1 = Point(-100, -100)
CP2 = Point(100, -100)
P2  = Point(200, 0)
@png begin
    circle.([P1, CP1, CP2, P2], 5, :fill)
    label.(["P1", "CP1", "CP2", "P2"], :N, [P1, CP1, CP2, P2], offset=10)
end 800 300 "images/bezier/fourpoints.png"

four points

The first and last points, P1 and P2, are the start and end of the line. The second point, CP1, controls the direction of the line as it leaves the first point, and the third point, CP2, determines the way the line approaches the fourth point. PostScript guru Don Lancaster (see footnote) uses the terms ‘influence point’ and ‘enthusiasm’—so the second influence point determines the enthusiasm with which the curve travels towards its final destination.

The graphics primitive curve() function draws a Bézier curve between P1 and P2, taking into account the positions of control points CP1 and CP2. It takes just three points, using the current position as the starting point.

using Colors
@png begin
    sethue("red")
    label.(["P1", "CP1", "CP2", "P2"], :N, [P1, CP1, CP2, P2])
    line(P1, CP1,:stroke)
    line(CP2, P2,:stroke)
    sethue("black")
    move(P1)
    curve(CP1, CP2, P2)
    strokepath()
    circle.([P1, CP1, CP2, P2], 3, :fill)
end 800 300 "images/bezier/curve.png"

simple bezier curve

The control points are like handles, controlling the shape of the curve. If you’ve used Adobe Illustrator or some other vector graphics software, you’ll be familiar with the idea of interactively dragging the handles around to get interesting curves. In Luxor, though, we sacrifice interactivity in favour of ruthless machine-driven automation. In the following animation, the control handles explore the geometry of a couple of hypotrochoids, while the helpless Bézier curve is pinned between them and forced into sinuous contortions:

animation gif

Can I make this in Adobe Illustrator? Hold my beer…

While you’re waiting, have a look at another animation; this is my artist’s impression of the De Casteljau algorithm dividing the control polygons around a Bézier curve as the parameter n moves from 0 to 1. The idea is that as p1 divides A to A1, p2 divides A1 to B1 and p3 divides B1 to B. So, pp1 divides p1 to p2, and pp2 divides p2 to p3. And you keep doing this until you can’t divide any more, and eventually the point P plots the course of the final Bézier curve. The red and blue parts of the curve show that this technique is also a good way to split a single Bézier curve into two separate ones, and the red and blue parts are separate control polygons.

De Casteljau

In Luxor, a BezierPathSegment type contains four 2D points, stored in the fields p1, cp1, cp2, and p2, and a BezierPath is an array of one or more of these BezierPathSegments. The drawbezierpath() function draws a BezierPath or BezierPathSegment, with similar results to the curve() function:

@png begin
    P1  = Point(-200, 0)
    CP1 = Point(-100, -100)
    CP2 = Point(100, -100)
    P2  = Point(200, 0)

    sethue("red")
    label.(["P1", "CP1", "CP2", "P2"], :N, [P1, CP1, CP2, P2])
    line(P1, CP1,:stroke)
    line(P2, CP2,:stroke)
    sethue("black")

    beziersegment = BezierPathSegment(P1, CP1, CP2, P2)
    circle.(beziersegment, 3, :fill)
    drawbezierpath(beziersegment, :stroke)
end 800 250 "images/bezier/drawbezierpath.png"

draw a Bezier path

This last function, drawbezierpath() is a typical Luxor drawing function, in that you can provide :fill as an alternative action to :stroke.

An easy way to make a BezierPath is to use makebezierpath() and supply a polygon. For example, let’s make a triangle with ngon() and use it as the skeleton for a new Bézier path:

@png begin
    sethue("red")
    triangle = ngon(O, 120, 3, vertices=true)
    poly(triangle, :stroke, close=true)

    sethue("black")
    bezierpath = makebezierpath(triangle)
    drawbezierpath(bezierpath, :stroke)
    circle.(triangle, 5, :fill)

    for bps in bezierpath
        circle.(bps, 3, :fill)
        line(bps[1], bps[2], :stroke)
        line(bps[3], bps[4], :stroke)
    end

    label.(["p1", "p2", "p3"], [:S, :NW, :E], triangle, offset=10)
end 800 400 "images/bezier/makebezierpath.png"

make bezier path

Here makebezierpath() converted the three points into an array of three separate Bézier path segments. The control points are positioned so that the curve flows freely from one segment to the segment.

Going straight

Bézier curves can have straight bits too. This animation shows the control points moving towards the points they’re controlling. When they merge, the Bézier path appears to become a series of straight lines:

function frame(scene, framenumber)
    background("white")
    eased_n = scene.easingfunction(framenumber, 0, 1, scene.framerange.stop)
    sethue("red")
    triangle = ngon(O, 120, 3, vertices=true)
    poly(triangle, :stroke, close=true)
    sethue("black")
    bezierpath = makebezierpath(triangle)
    newbezierpath = BezierPath()
    for bps in bezierpath
        push!(newbezierpath,
            BezierPathSegment(bps.p1,
             between(bps.cp1, bps.p1, eased_n),
             between(bps.cp2, bps.p2, eased_n),
             bps.p2))
    end
    for bps in newbezierpath
        circle.(bps, 3, :fill)
        line(bps.p1, bps.cp1, :stroke)
        line(bps.cp2, bps.p2, :stroke)
    end
    drawbezierpath(newbezierpath, :stroke)
end

width, height = (256, 256)

bedtimemovie = Movie(width, height, "bedtimemovie")
# probably requires ffmpeg installed
animate(bedtimemovie, [Scene(bedtimemovie, frame, 1:150, easingfunction=easeoutquad)], creategif=true, framerate=30, pathname="images/bezier/bedtimemovie.gif")

bedtime for Béziers

This is a process known to very young Adobe Illustrator users as ‘putting the Bézier handles to bed’.

It’s also fun to move the control points somewhere else. Here, they’re multiplied by 2, while the first and last points are left unchanged:

@png begin
    sethue("red")
    triangle = ngon(O, 100, 3, vertices=true)
    poly(triangle, :stroke, close=true)

    sethue("purple")
    bezierpath = makebezierpath(triangle)
    drawbezierpath(bezierpath, :stroke)
    sethue("magenta")
    for bps in bezierpath
        bps = BezierPathSegment(bps[1], 2bps[2], 2bps[3], bps[4])
        drawbezierpath(bps, :stroke)
    end

end 800 400 "images/bezier/movecontrolpoints.png"

moving control points

We could make even more copies, multiplying the control points each time through:

@png begin
    setline(.5)
    pgon = ngon(O, 100, 3, vertices=true)
    bezierpath = makebezierpath(pgon)

    for bpseg in bezierpath
        for i in 1:20
            bpseg.cp1 *= 1.05
            bpseg.cp2 *= 1.05
            drawbezierpath(bpseg, :stroke)
        end
    end
end 800 400 "images/bezier/movecontrolpointsmore.png"

moving multiple control points

Try changing the initial triangle to a pentagon or heptagon by changing the 3 in ngon(). And try multiplying by less than 1.05 too…

Should you ever want to draw Bézier curves “the hard way”, try this. The standard Bézier function is available in Luxor as bezier(), and we could draw a small hue-varying circle at each point:

@png begin
    setopacity(0.6)
    P1  = Point(-250, 0)
    CP1 = Point(-200, -150)
    CP2 = Point(200, -150)
    P2  = Point(250, 0)
    circle.([P1, CP1, CP2, P2], 3, :fill)

    for u in 0.0:0.025:1.0
        sethue(Colors.HSV(rescale(u, 0, 1, 0, 360), 1, 1))
        pt = bezier(u, P1, CP1, CP2, P2)
        circle(pt, 5, :fill)
    end

end 800 400 "images/bezier/thehardway.png"

bezier the hard way

What happens if you step u from say, -25.0 to 25.0? And is that supposed to happen? (Spoiler: the curve doubles back and shoots off to (±∞/∞).)

Out of control

Suppose you wanted to draw a Bézier curve but you didn’t know where the control points were, but you did know a couple of points that should lie on the line? Again, Luxor has the answer for you, in the form of a function called bezierfrompoints(). You supply four points, and the function returns the points of the Bézier curve that passes through them.

@png begin
    sethue("grey80")
    box(O, 120, 100, :stroke)
    corners = box(O, 120, 100, vertices=true)
    label.(["1", "2","3", "4"], :se, corners, offset=10)
    sethue("red")
    circle.(corners, 5, :fill)

    startpoint, cp1, cp2, endpoint  = bezierfrompoints(corners...)

    sethue("purple")
    line(startpoint, cp1, :stroke)
    line(endpoint, cp2, :stroke)

    sethue("turquoise")
    circle.([cp1, cp2, startpoint, endpoint], 3, :fill)

    sethue("blue")
    drawbezierpath(BezierPathSegment(startpoint, cp1, cp2, endpoint), :stroke)
end 800 300 "images/bezier/bezierfrompoints.png"

bezier through four points

Note that this function returns all four points, but of course you already knew the first and last ones (in corners), you just wanted the two control points.

On the right path?

In Luxor the three main ways to make graphic lines and shapes are: paths, polygons, and BezierPaths. Paths are the fundamental building blocks of graphics, consisting of one or more sequences of straight and Bézier curves. A polygon is an array of points, which will be converted to a path when you draw it. And BezierPaths consist of a list of BezierPathSegments, which will also be converted to ordinary paths when they’re drawn.

It’s useful to be able to convert between the different types:

Function Converts
makebezierpath(pgon) polygon to BezierPath
pathtopoly() current path to array of polygons
pathtobezierpaths() current path to array of BezierPaths
beziertopoly(bpseg) BezierPathSegment to polygon
bezierpathtopoly(bezierpath) BezierPath to polygon
bezierfrompoints(p1, p2, p3, p4) convert four points to BezierPathSegment

Can’t draw a circle?

Having witnessed the ability of a simple Bézier curve to adopt so many shapes, it’s a bit surprising to find out that you can’t use a Bézier curve to draw a circle. Well, it can do a very good impression of one, of course, but mathematically it can’t produce a purely circular curve. You can see this if you draw a very large circle and a very large matching Bézier segment. For a circle with radius 10 meters (and that’s one big PDF), the discrepancy between the pure circle and the Bézier approximation isn’t much bigger than the size of this period/full stop.←

We non-scientists are lucky in not having to worry about errors of this magnitude…

a 10 metre circle

In the above picture, the red circle is made by circle(), the green one is made with Bézier curves (as used by circlepath() and ellipse()).

“0.55228474983 is the magic number”

To draw approximate circles using Bézier curves, you need to know the magic number, usually called kappa, which has the value 0.552284…. We’re looking at a Bézier curve pretending to be a circular quadrant, with control points positioned at a certain distance kappa from the end points. But how far? What value of kappa—what length of handle—will give us the most circular curve?

A picture of the problem, using a deliberately not-very-good guess at a value for kappa of 0.5:

@png begin
    scale(3)
    fontsize(7)
    translate(0, -50)
    setline(1)

    radius = 100

    kappa_guess = .5

    sethue("red")
    arc(O, radius, 0, pi/2, :stroke)

    sethue("black")
    quadrant = BezierPathSegment(
        Point(0, radius),
        Point(kappa_guess * radius, radius),
        Point(radius, kappa_guess * radius),
        Point(radius, 0))
    circle.(quadrant, 1, :fill)
    arrow(quadrant.p1, quadrant.cp1, arrowheadlength=5)
    arrow(quadrant.p2, quadrant.cp2, arrowheadlength=5)

    drawbezierpath(quadrant, :stroke)

    arrow(O, bezier(0.5, quadrant...))

    label.(["(0, 1)", "(k, 1)", "(1, k)", "(1, 0)"], [:S, :S, :E, :E], quadrant)
    label("kappa_guess = 0.5", :S, midpoint(quadrant.p1, quadrant.cp1), offset=15)
end 800 450 "images/bezier/kappaguess.png"

guessing at kappa

To solve this, we’ll define a function that takes a value for kappa and works out the radius at the center of the Bézier curve.

_radiusofbezier(kappa) = norm(O,
    bezier(0.5, BezierPathSegment(
        Point(0, 1),
        Point(kappa * 1, 1),
        Point(1, kappa * 1),
        Point(1, 0)
        )...))

radiusofbezier(x) = _radiusofbezier(x) - 1

And we’ll use one of the many excellent packages in JuliaMath, Roots.jl, to find the zero point:

using Roots
find_zero(radiusofbezier, (0, 1), verbose=true)

Results of univariate zero finding:
* Converged to: 0.5522847498307936
* Algorithm: Roots.Bisection64()
* iterations: 59
* function evaluations: 61
* stopped as |f(x_n)|  max(δ, max(1,|x|)⋅ϵ) using δ = atol, ϵ = rtol
Trace:
x_0 =  0.0000000000000000,   f(x_0) = -0.2928932188134524
x_1 =  0.0000000000000000,   f(x_1) = -0.2928932188134524
x_2 =  0.0000000000000000,   f(x_2) = -0.2928932188134524
x_3 =  0.0000000000000000,   f(x_3) = -0.2928932188134524
x_4 =  0.0000000000000000,   f(x_4) = -0.2928932188134524
x_5 =  0.0000000002401066,   f(x_5) = -0.2928932186861167
x_6 =  0.0000154972076416,   f(x_6) = -0.2928850001779929
x_7 =  0.0039367675781250,   f(x_7) = -0.2908054325256169
x_8 =  0.0627441406250000,   f(x_8) = -0.2596181133267076
x_9 =  0.2504882812500000,   f(x_9) = -0.1600517471037238
x_10 =  0.5004882812500000,  f(x_10) = -0.0274692256312462
x_11 =  0.5004882812500000,  f(x_11) = -0.0274692256312462
x_12 =  0.5004882812500000,  f(x_12) = -0.0274692256312462
x_13 =  0.5004882812500000,  f(x_13) = -0.0274692256312462
x_14 =  0.5317077636718750,  f(x_14) = -0.0109125948370147
x_15 =  0.5473175048828125,  f(x_15) = -0.0026342794398989
x_16 =  0.5473175048828125,  f(x_16) = -0.0026342794398989
x_17 =  0.5512199401855469,  f(x_17) = -0.0005647005906200
x_18 =  0.5512199401855469,  f(x_18) = -0.0005647005906200
x_19 =  0.5521955490112305,  f(x_19) = -0.0000473058783003
x_20 =  0.5521955490112305,  f(x_20) = -0.0000473058783003
x_21 =  0.5521955490112305,  f(x_21) = -0.0000473058783003
x_22 =  0.5521955490112305,  f(x_22) = -0.0000473058783003
x_23 =  0.5522565245628357,  f(x_23) = -0.0000149687087803
x_24 =  0.5522565245628357,  f(x_24) = -0.0000149687087803
x_25 =  0.5522717684507370,  f(x_25) = -0.0000068844164003
x_26 =  0.5522793903946877,  f(x_26) = -0.0000028422702103
x_27 =  0.5522832013666630,  f(x_27) = -0.0000008211971153
x_28 =  0.5522832013666630,  f(x_28) = -0.0000008211971153
x_29 =  0.5522841541096568,  f(x_29) = -0.0000003159288415
x_30 =  0.5522846304811537,  f(x_30) = -0.0000000632947047
x_31 =  0.5522846304811537,  f(x_31) = -0.0000000632947047
x_32 =  0.5522847495740280,  f(x_32) = -0.0000000001361705
x_33 =  0.5522847495740280,  f(x_33) = -0.0000000001361705
x_34 =  0.5522847495740280,  f(x_34) = -0.0000000001361705
x_35 =  0.5522847495740280,  f(x_35) = -0.0000000001361705
x_36 =  0.5522847495740280,  f(x_36) = -0.0000000001361705
x_37 =  0.5522847495740280,  f(x_37) = -0.0000000001361705
x_38 =  0.5522847495740280,  f(x_38) = -0.0000000001361705
x_39 =  0.5522847495740280,  f(x_39) = -0.0000000001361705
x_40 =  0.5522847495740280,  f(x_40) = -0.0000000001361705
x_41 =  0.5522847498066312,  f(x_41) = -0.0000000000128140
x_42 =  0.5522847498066312,  f(x_42) = -0.0000000000128140
x_43 =  0.5522847498066312,  f(x_43) = -0.0000000000128140
x_44 =  0.5522847498066312,  f(x_44) = -0.0000000000128140
x_45 =  0.5522847498211689,  f(x_45) = -0.0000000000051041
x_46 =  0.5522847498284378,  f(x_46) = -0.0000000000012493
x_47 =  0.5522847498284378,  f(x_47) = -0.0000000000012493
x_48 =  0.5522847498302550,  f(x_48) = -0.0000000000002855
x_49 =  0.5522847498302550,  f(x_49) = -0.0000000000002855
x_50 =  0.5522847498307093,  f(x_50) = -0.0000000000000447
x_51 =  0.5522847498307093,  f(x_51) = -0.0000000000000447
x_52 =  0.5522847498307093,  f(x_52) = -0.0000000000000447
x_53 =  0.5522847498307660,  f(x_53) = -0.0000000000000145
x_54 =  0.5522847498307660,  f(x_54) = -0.0000000000000145
x_55 =  0.5522847498307802,  f(x_55) = -0.0000000000000070
x_56 =  0.5522847498307873,  f(x_56) = -0.0000000000000032
x_57 =  0.5522847498307909,  f(x_57) = -0.0000000000000013
x_58 =  0.5522847498307927,  f(x_58) = -0.0000000000000004
x_59 =  0.5522847498307936,  f(x_59) =  0.0000000000000000
0.5522847498307936

So that’s kappa.

Alternatively, we can use algebraic sorcery and the parametric equation for the Bézier cubic function to conjure the value from the following incantations:

I suspect most applications simply hard-code the magic number 0.552284… directly.

The following code shows how we could animate the process of changing the length of the handles by changing the value of kappa. Blink and you’ll miss the sweet spot, though, and the movement of the handles is imperceptible, so perhaps this isn’t the best way of illustrating the construction.

function frame(scene, framenumber)
    background("white")
    eased_n = scene.easingfunction(framenumber, 0, 1, scene.framerange.stop)

    radius = 250
    setline(1)
    sethue("black")
    kappa = rescale(eased_n, 0, 1, 0.525, 0.575)

    handlelength = (radius/2)  * kappa

    left   = Point(-radius/2, 0)
    top    = Point(0,        -radius/2)
    right  = Point(radius/2,  0)
    bottom = Point(0,         radius/2)

    line(left, left + (0, -handlelength), :stroke)
    line(top + (-handlelength, 0), top, :stroke)
    circle.([left, top, right, bottom], 2, :fill)
    circle.([left + (0, -handlelength), top + (-handlelength, 0)], 2, :fill)

    line(left, left + (0, -handlelength), :stroke)
    line(top + (-handlelength, 0), top, :stroke)
    circle.([left, top, right, bottom], 2, :fill)
    circle.([left + (0, -handlelength), top + (-handlelength, 0)], 2, :fill)

    move(left)
    curve(left   + (0, -handlelength),
          top    + (-handlelength, 0),
          top)
    curve(top    + (handlelength, 0),
          right  + (0, -handlelength),
          right)
    curve(right  + (0, handlelength),
          bottom + (handlelength, 0),
          bottom)
    curve(bottom + (-handlelength, 0),
          left   + (0, handlelength),
          left)

    if isapprox(kappa, 0.55228, atol=0.001)
        sethue("red")
        fillpath()
        sethue("black")
    else
        strokepath()
    end
    fontface("Menlo")
    fontsize(20)
    text(string(round(kappa, 4)), halign=:center)
end
width, height = (400, 400)
circlefrombezier = Movie(width, height, "circlefrombezier")
animate(circlefrombezier, [Scene(circlefrombezier, frame, 1:150, easingfunction=easeoutquad)], creategif=true, framerate=30, pathname="images/bezier/circlefrombezier.gif")

kappa value animation

Luxor already provides a Luxor.circlepath() function that uses four Bézier curve segments to build a path that draws a circle. The main advantage of using this instead of the default arc-based circle is that it’s easier to build circles with holes:

@png begin
    background("moccasin")
    sethue("darkorange")
    circlepath(O, 200, :path)

    newsubpath()
    circlepath(O + (100, -40), 70, :path, reversepath=true)

    newsubpath()
    circlepath(O - (100, 40), 70, :path, reversepath=true)

    fillpreserve()
    setline(8)
    sethue("skyblue")
    strokepath()
end 800 500 "images/bezier/circularpaths.png"

circular paths

Detour into Typomania

In practice, not being able to draw perfect circles with Bezier curves isn’t a big problem. Font designers (you knew this was going to go ‘all typographical’ sooner or later) are big users of Bézier curves, and typically they don’t like drawing perfect circles anyway, because of the optical corrections required to make things ‘look right’.

There are a number of optical illusions that demonstrate that the human eye and the brain don’t always see reality accurately. The following is one of the simplest, but most people would be prepared to bet that the horizontal bar is thicker and shorter than the vertical bar. I had to draw a grid to double-check…

@png begin
    background("white")

    sethue("gray90")
    g = GridRect(BoundingBox()[1], 15, 15, 30 * 20, 30 * 20)
    for i in 1:1600
        circle(nextgridpoint(g), 0.5, :fill)
    end

    b = box(O - (15, -60), O + (15, -60), vertices=true)

    sethue("black")
    poly(b, :fill)
    translate(0, -90)
    rotate(pi/2)
    poly(b, :fill)

end 300 300 "images/bezier/opticalillusion.png"

simple illusion

It’s something to do with how our eyes, set side by side and trained to move side to side horizontally with great speed and precision, underestimate width. Type designers spend much of their time adjusting the relative widths and thicknesses of letter shapes so that illusions like this are compensated for in advance. What happens if you turn your display on its side (apart from possibly spilling your coffee)?

Here’s a short script to examine the bounding boxes of the inner and outer loops of the letter ‘o’ in various fonts.

function drawletteraspect(pos, font, thefontsize, str)
    @layer begin
        translate(pos)
        translate(-thefontsize/3, 0)
        fontsize(thefontsize)
        fontface(font)
        sethue("gray40")
        setline(0.5)
        text(str)
        newpath()
        textpath(str)
        o = pathtopoly()
        @layer begin
            fontsize(thefontsize/20)
            text(font, Point(0, thefontsize/15))
            sethue("magenta")
            fontface("Cousine-Bold")
            n = 0
            for i in o
                b = BoundingBox(i)
                text(string(round(boxaspectratio(b), 4)), midpoint(b...) + (0, n), halign=:center)
                n += thefontsize/20
                box(b, :stroke)
            end
        end
    end
end

@png begin
    fonts = [
        "ArialUnicodeMS",
        "AvantGarde-Book",
        "BauhausLT-Light",
        "CircularStd-Book",
        "Cousine",
        "FuturaLT-Book",
        "GillSansMTPro-Book",
        "GothamBlack",
        "Helvetica",
        "KabelLTStd-Light",
        "ProximaNova-Bold",
        "Times-Roman",
    ]
    tiles = Tiler(800, 400, 3, 4, margin=40)
    for (pos, n) in tiles
        drawletteraspect(pos, fonts[mod1(n, length(fonts))], 150, "o")
    end
end 800 400 "images/bezier/drawletteraspect.png"

aspect ratios of letter O

The fonts on your system will be different, of course.

The story of “o”

I wondered which fonts used the most circular circles for the letter “o”. This script wanders through all the fonts registered with Fontconfig and stores the bounding boxes of the letter “o”, then finds “the most circular”.

using Fontconfig, Luxor, DataFrames

function buildfontlist()
    fonts = []
    for font in Fontconfig.list()
        families = Fontconfig.format(font, "%{family}")
        for family in split(families, ",")
            push!(fonts, family)
        end
    end
    filter!(font -> !startswith(font, "."), fonts)
    filter!(font -> !contains(font, r"Wing.*"), fonts)
    filter!(font -> !contains(font, r".*Extra"), fonts)
    filter!(font -> !contains(font, r".*Expert"), fonts)
    return sort(unique(fonts))
end

mutable struct Score
    fontname::String
    bbox1::BoundingBox
    bbox2::BoundingBox
end

scores = Score[]

# we have to have a drawing open, otherwise textpath() doesn't know where to turn

Drawing(100, 100, "/tmp/emptydrawing.png")
fontlist = unique(buildfontlist())
fontsize(20)

for font in fontlist
    fontface(font)
    newpath() # !!!!
    textpath("o")
    o = pathtopoly()
    bboxes = BoundingBox[]
    for i in o
        push!(bboxes, BoundingBox(i))
    end
    if length(bboxes) == 2
        push!(scores, Score(font, bboxes[1], bboxes[2]))
    end
end
finish() # unused drawing

To analyse the results, we’ll put everything into a DataFrame.

df = DataFrame(
    Fontname = String[],  # font name
    Outer    = Float64[], # outer aspect ratio
    Inner    = Float64[], # inner loop aspect ratio
    Total    = Float64[],
    Mean     = Float64[]
    )

for f in scores
    push!(df, [f.fontname,
               boxaspectratio(f.bbox1),
               boxaspectratio(f.bbox2),
               1 - (boxaspectratio(f.bbox1) + boxaspectratio(f.bbox2))/2,
               mean([boxaspectratio(f.bbox1),
                     boxaspectratio(f.bbox2)])
              ])
end

sort!(df, :Total, lt = (a, b) -> abs(a) < abs(b))

1018×5 DataFrames.DataFrame
 Row   Fontname                         Outer     Inner    Total         Mean     
├──────┼─────────────────────────────────┼──────────┼─────────┼──────────────┼──────────┤
 1     Kabel LT Std Light               1.0       1.0      0.0           1.0      
 2     ITC Avant Garde Std XLt          0.9965    1.0035   -2.04211e-7   1.0      
 3     Kabel LT Std Black               1.0       1.00101  -0.000503525  1.0005   
 4     Gill Sans                        0.972742  1.02949  -0.00111546   1.00112  
 5     Gill Sans MT Pro Medium          0.972742  1.02949  -0.00111546   1.00112  
 6     ITC Avant Garde Gothic Std       0.993158  1.01428  -0.00371842   1.00372  
 7     ITC Avant Garde Std Bk           0.993158  1.01428  -0.00371842   1.00372  
 8     ITC Avant Garde Std Md           0.991246  1.021    -0.00612158   1.00612  
 9     Avenir Heavy                     0.965208  1.02078  0.00700837    0.992992 

 1009  Akzidenz-Grotesk BQ Condensed A  2.30537   5.19277  -2.74907      3.74907  
 1010  AmplitudeComp-Ultra              1.52721   6.41949  -2.97335      3.97335  
 1011  Industria LT Std Solid           2.62973   5.72823  -3.17898      4.17898  
 1012  Bernard MT Condensed             1.5877    7.55733  -3.57252      4.57252  
 1013  Impact                           1.4851    7.6678   -3.57645      4.57645  
 1014  HeadLineA                        1.6       8.34783  -3.97391      4.97391  
 1015  Univers LT Std 39 Thin UltraCn   3.69307   6.81657  -4.25482      5.25482  
 1016  Helvetica LT Std ExtCompressed   1.9595    8.82251  -4.39101      5.39101  
 1017  Gill Sans MT Pro Bold ExtCond    2.70499   9.1913   -4.94815      5.94815  
 1018  Helvetica LT Std UltCompressed   2.46103   9.71489  -5.08796      6.08796  

Few typefaces are perfectly circular. Even Circular, LineTo’s trendy geometric sans typeface, isn’t perfectly circular.

Most of the least circular ones, according to this rough examination, are in the Condensed and Compressed sections of the font libraries. No surprise, Sherlock.

display(df[end-10:end, :])

11×5 DataFrames.DataFrame
 Row  Fontname                         Outer    Inner    Total     Mean    
├─────┼─────────────────────────────────┼─────────┼─────────┼──────────┼─────────┤
 1    Helvetica LT Std Compressed      5.72485  1.53453  -2.62969  3.62969 
 2    Akzidenz-Grotesk BQ Condensed A  2.30537  5.19277  -2.74907  3.74907 
 3    AmplitudeComp-Ultra              1.52721  6.41949  -2.97335  3.97335 
 4    Industria LT Std Solid           2.62973  5.72823  -3.17898  4.17898 
 5    Bernard MT Condensed             1.5877   7.55733  -3.57252  4.57252 
 6    Impact                           1.4851   7.6678   -3.57645  4.57645 
 7    HeadLineA                        1.6      8.34783  -3.97391  4.97391 
 8    Univers LT Std 39 Thin UltraCn   3.69307  6.81657  -4.25482  5.25482 
 9    Helvetica LT Std ExtCompressed   1.9595   8.82251  -4.39101  5.39101 
 10   Gill Sans MT Pro Bold ExtCond    2.70499  9.1913   -4.94815  5.94815 
 11   Helvetica LT Std UltCompressed   2.46103  9.71489  -5.08796  6.08796 

t = Table(2, 5, 170, 140)
@png begin
    for (n, f) in enumerate(df[end-9:end, 1])
        fontsize(60)
        fontface(f)
        text("O", t[n], halign=:center)
        fontsize(14)
        text(f, t[n] + (0, 40), halign=:center)
    end
end 800 250 "images/bezier/leastcircularfonts.png"

some of the least circular fonts

The most circular ones on my computer are the classic geometric sans serif fonts.

display(df[1:10, :])

10×5 DataFrames.DataFrame
 Row  Fontname                    Outer     Inner    Total         Mean     
├─────┼────────────────────────────┼──────────┼─────────┼──────────────┼──────────┤
 1    Kabel LT Std Light          1.0       1.0      0.0           1.0      
 2    ITC Avant Garde Std XLt     0.9965    1.0035   -2.04211e-7   1.0      
 3    Kabel LT Std Black          1.0       1.00101  -0.000503525  1.0005   
 4    Gill Sans                   0.972742  1.02949  -0.00111546   1.00112  
 5    Gill Sans MT Pro Medium     0.972742  1.02949  -0.00111546   1.00112  
 6    ITC Avant Garde Gothic Std  0.993158  1.01428  -0.00371842   1.00372  
 7    ITC Avant Garde Std Bk      0.993158  1.01428  -0.00371842   1.00372  
 8    ITC Avant Garde Std Md      0.991246  1.021    -0.00612158   1.00612  
 9    Avenir Heavy                0.965208  1.02078  0.00700837    0.992992 
 10   DIN 30640 Std               0.966066  1.01366  0.010139      0.989861 

t = Table(2, 5, 160, 140)
@png begin
    for (n, f) in enumerate(df[1:10, 1])
        fontsize(50)
        fontface(f)
        text("O", t[n], halign=:center)
        fontsize(11)
        text(f, t[n] + (0, 50), halign=:center)
    end
end 800 300 "images/bezier/mostcircularfonts.png"

the most circular fonts

In the ‘most circular o’ fonts, there are a few examples from Rudolf Koch’s Kabel family.

kabel = df[[contains(fontname, r".*Kabel.*") for fontname in df[:Fontname]], :]

9×5 DataFrames.DataFrame
 Row  Fontname              Outer     Inner    Total         Mean     
├─────┼──────────────────────┼──────────┼─────────┼──────────────┼──────────┤
 1    Kabel LT Std Light    1.0       1.0      0.0           1.0      
 2    Kabel LT Std Black    1.0       1.00101  -0.000503525  1.0005   
 3    ITC Kabel Std Medium  0.976604  0.99728  0.0130582     0.986942 
 4    Kabel LT Std          1.02945   1.0      -0.014727     1.01473  
 5    Kabel LT Std Book     1.02945   1.0      -0.014727     1.01473  
 6    ITC Kabel Std Ultra   0.961288  1.07731  -0.0193006    1.0193   
 7    ITC Kabel Std         1.02357   1.02993  -0.0267482    1.02675  
 8    ITC Kabel Std Book    1.02357   1.02993  -0.0267482    1.02675  
 9    ITCKabel              1.02357   1.02993  -0.0267482    1.02675  

# that's v0.6 syntax... v0.7 is df[[occursin(r".*Kabel.*", fname) for fname  in df[:Name]], :]

t = Table(2, 5, 160, 140)
@png begin
    for (n, f) in enumerate(kabel[:, 1])
        fontsize(50)
        fontface(f)
        text("O", t[n], halign=:center)
        fontsize(11)
        text(f, t[n] + (0, 50), halign=:center)
    end
end 800 300 "images/bezier/kabelfonts.png"

Kabel fonts

Here’s Wikipedia:

Kabel belongs to the “geometric” style of sans-serifs, which was becoming popular in Germany at the time of Kabel’s creation [1920s]. Based loosely on the structure of the circle and straight lines, it nonetheless applies a number of unusual design decisions, such as a delicately low x-height (although larger in the bold weight), a quirky tilted ‘e’ and irregularly angled terminals, to add delicacy and an irregularity suggesting stylish calligraphy, of which Koch was an expert.

Eye magazine isn’t convinced by the apparent geometrical precision of Kabel:

its eccentricities reveal Koch’s unwavering expressionistic and humanist instincts

But at least the ‘o’s are circular…

Curvy

Moving back to Bézier curves, you can summarize the behaviour of a Bézier curve at a point by finding the curvature. The Luxor function beziercurvature() returns a number, called kappa, that varies and flips from positive to negative as the Bézier path varies in ‘curviness’. We’re using the first and second derivatives to find the kappa value:

The following drawcurvature() function uses the kappa value to work out the slope, and draw a perpendicular to the curve, with lengths varying according to the value of kappa, indicating the way the curvature changes.

(This is another kappa, by the way, no relation to the Bézier circularity kappa. Or indeed to the Lancia Kappa…

Lancia Kappa from Wikimedia Commons

Photograph by Rudolf Stricker

…or to many of the other things called kappa, such as the curvature of the universe, the torsional constant of an oscillator, Einstein’s constant of gravitation, the coupling coefficient in magnetostatics—and that’s just in physics.)

function drawcurvature(bezpathsegment::BezierPathSegment;
        scale=500,
        stepby=0.01,
        colors=["purple", "green"])
    p1, cp1, cp2, p2 = bezpathsegment
    for u in 0:stepby:1
        pt = bezier(u, p1, cp1, cp2, p2)
        kappa = beziercurvature(u, p1, cp1, cp2, p2)
        # comb lines are perpendiculars
        @layer begin
            pt1 = bezier′(u, p1, cp1, cp2, p2)
            # slope of curve at t = y′(t)/x′(t)
            # normal is 1/slope
            s = -1/(pt1.y/pt1.x)
            sign(kappa) >= 0 ? sethue(colors[1]) : sethue(colors[2])
            len = scale * abs(kappa)
            line(pt - polar(len, atan(s)), pt + polar(len, atan(s)), :stroke)
        end
    end
end

Here it is in action:

@png begin
    p1  = Point(-150, 0)
    cp1 = Point(-50, -200)
    cp2 = Point(150, 200)
    p2  = Point(250, -150)
    setline(0.1)
    for u in 0:0.025:1
        pt = bezier(u, p1, cp1, cp2, p2)
        circle(pt, 2, :fill)
    end
    sethue("red")
    circle.([p1, cp1, cp2, p2], 3, :fill)
    setline(1)
    line(p1, cp1, :stroke)
    line(p2, cp2, :stroke)
    sethue("black")
    fontsize(6)
    setline(0.5)
    drawcurvature(BezierPathSegment(p1, cp1, cp2, p2),
        scale=1500,
        stepby=0.005)
end 600 500 "images/bezier/drawcurvaturecomb.png"

drawing the curvature comb

Users of CAD systems (particularly industrial designers) like to display these curvature combs (in both 2D and 3D) to make sure that their shapes don’t introduce noticeably abrupt (and possible weak) transitions. Type designers use them too, but as usual they like to let their eyes have the final say.

@png begin
    background("gray30")
    fontface("Times-BoldItalic")
    translate(-250, 80)
    fontsize(300)
    textpath("julia")
    bpaths = pathtobezierpaths()
    sethue("gray60")
    for bpath in bpaths
        for bpseg in bpath
            setline(0.5)
            drawbezierpath(bpseg, :stroke)
            setline(0.25)
            drawcurvature(bpseg,
                scale=180,
                stepby=0.05,
                colors=["gold", "cyan"])
        end
    end
end 800 500 "images/bezier/typecomb.png"

curvature comb for the g

The dots over the ‘i’ and ‘j’ don’t look like perfect circles; these are defined by 8 points, not 4, for some reason.

Osculate my Béziers

The value of kappa is typically very small, so the radius of curvature, which is defined as , can become very large. This radius value defines a circle that just touches the curve and follows the curvature at that point. Mathematicians, in typically romantic mood, call it the osculating circle, osculate being from the Latin noun osculum, meaning “kiss”.

Drawing osculating circles can be a challenge; they grow very large when the curve looks flat. To be honest it’s a tricky diagram to style up, and there’s a lot of information that it would be cool to keep and a shame to throw away. There’s quite a bit of osculation going on here…

@png begin
    p1  = Point(-400, 60)
    cp1 = Point(-300, -400)
    cp2 = Point(250, 300)
    p2  = Point(200, -100)
    background("pink")
    scale(0.75)
    for u in 0:0.001:1
        pt = bezier(u, p1, cp1, cp2, p2)
        circle(pt, 4, :fill)
    end
    sethue("red")
    circle.([p1, cp1, cp2, p2], 5, :fill)
    setline(2)
    line(p1, cp1, :stroke)
    line(p2, cp2, :stroke)
    setopacity(0.1) ; setline(1.5) ; fontsize(5)
    for pointofcontact in 0.0:0.015:1.0
        kappa = beziercurvature(pointofcontact, p1, cp1, cp2, p2)
        pt = bezier(pointofcontact, p1, cp1, cp2, p2)
        @layer begin
            setopacity(1.0)
            sethue("black")
            circle(pt, 5, :fill)
        end
        pt1 = bezier′(pointofcontact, p1, cp1, cp2, p2)
        s = -1/(pt1.y/pt1.x) # normal is 1/slope
        radius = abs(1/kappa)
        centrepoint = pt + polar(radius, atan(s))
        if radius < 600
            sign(kappa) >= 0 ? sethue("plum1") : sethue("mediumpurple")
            circle(centrepoint, radius, :fill)
            @layer begin
                setopacity(1)
                setline(0.25)
                sethue("black")
                fontsize(12)
                circle(centrepoint, radius, :stroke)
                circle(centrepoint, 2, :fill)
                line(centrepoint, pt, :stroke)
                text(string(round(1/kappa, 2)), centrepoint)
            end
        else
            text("too big", pt)
        end
    end
end 800 800 "images/bezier/osculate.png"

osculating

A blot on the landscape

Bézier curves became popular because they allow us to specify gently curving shapes that would be impractical and cumbersome to specify with circular arcs. The Comprehensive Taxonomy of Irregular Amoeboid Shapes is perhaps still waiting to be written, but here’s my contribution to the ‘random ink blot’ chapter. Sometimes you’ll get lucky and get a nice one.

function blot(;pos=O, radius=50, npoints=10, action=:fill)
    center = pos
    pts = ngon(pos, radius, npoints, vertices=true)
    for i in 1:length(pts)
        pts[i] = pts[i] + (rand(-radius/3:radius/3), rand(-radius/3:radius/3))
    end
    bezpath = makebezierpath(pts)
    for n in 1:length(bezpath)
        bezseg = bezpath[mod1(n, length(bezpath))]
        if isodd(n)
            j = rand(1.4:0.1:7)
            bezseg.cp1 = between(center, bezseg.p1, j)
            bezseg.cp2 = between(center, bezseg.p2, j)
        else
            j = rand(0.25:0.1:0.7)
            bezseg.cp1 = between(center, bezseg.p1, j)
            bezseg.cp2 = between(center, bezseg.p2, j)
        end
    end
    drawbezierpath(bezpath, action)
    return bezpath
end

@png begin
    sethue("midnightblue")
    blot(npoints=10)
    circle.([Point(rand(-200:200), rand(-200:200)) for i in 1:15], rand(1:0.1:15, 15), :fill)
end 800 500 "images/bezier/blot.png"

blot

The control handles of alternate points are positioned on a line from the center to the point.

You could analyse the curvature of these blobs, if you really wanted to, using the drawcurvature() function from earlier:

@png begin
    pl = box(BoundingBox(), vertices=true)
    mesh1 = mesh(pl, ["purple", "blue", "yellow", "orange"])
    setmesh(mesh1)
    paint()
    setopacity(0.8)
    sethue("orange")
    setline(5)
    b = blot(radius=80, npoints=8, action=:stroke)
    drawbezierpath(b, :fill)
    setline(.25)
    drawcurvature.(b, colors=["white", "black"])
end 800 500 "images/bezier/blobcomb.png"

blob comb

A brush in the rough

Not all lines generated by computers have to be rigidly straight and precise. What if we could easily make graphics that are a bit more relaxed in style, rather than the rigid CAD-like (and yes, awesome) precision graphics we’re used to?

@png begin
    t = Table(1, 3)
    box(t[1], 160, 160, :stroke)
    boxv = box(t[3], 160, 160, vertices=true)
    foreach((f, t) ->
        brush(f, t, 3),
           boxv[1:4],
           boxv[mod1.(2:5, 4)])
end 800 300 "images/bezier/lines.png"

rough lines

“A line is a breadthless length.” (Euclid)

The idea here is that a single line between two points is replaced with some BezierPathSegments that together define a shape that can vary in thickness along its length. This shape can be then filled, and is independent of the set line thickness.

@png begin
    setline(60) # <— pointless
    brush(O - (250, 0), O + (250, 0), 4, strokes=1)
end 800 50  "images/bezier/brushline1.png"

brush line

The experimental brush() function is bristling with built-in randomness, so you never know what you’re going to get. There are some control knobs available which you can play with.

@png begin
    brush(O - (250, 0), O + (250, 0), 5,
        strokes=15,
        twist=4)
end 800 50  "images/bezier/brushline2.png"

brush line 2

To be honest, I think it’s a bit daft to abandon the machine-like precision that our graphics software usually gives us for this variable hand-made look.

using Iterators
@png begin
    background("gray20")
    sethue("ivory")
    for ss in subsets(ngon(O, 250, 13, vertices=true), 2)
        # line(ss[1], ss[2], :stroke) # use this line for accurate graphics
        brush(ss[1], ss[2], 4,
            strokes=1,
            minwidth=0.001,
            maxwidth=0.05,
            lowhandle=0.1,
            highhandle=.2,
            randomopacity=true)
    end
end 800 600  "images/bezier/linepartition.png"

line partition

Today, our relationship with mechanical production and product design is inconsistent; some of the things we desire we want to be hand-made, but others we’d prefer to be machine-made. Only the most expensive cars claim to be made ‘by hand’; the cheaper models flaunt their nanometre precision instead. Hipsters seek the authentic analogue roughness of the products of the Second Industrial age, but are secretly grateful for and rely on the smooth precisions afforded by the Third. Handmade shoes yes, handmode iPhones, no.

There are a number of plotting packages that offer a hand-drawn aesthetic—this site is web-based. So you can definitely announce your latest scientific discovery using XKCD-style presentation graphics. There are XKCD-styling kits for most of the software used by people who have heard of XKCD.

@png begin
    background("bisque")
    sethue("gray60")
    setopacity(0.1)
    blot(pos=O + (-10, -50), radius=100)
    sethue("springgreen2")
    brush(O + (20, 270), O + (0, -270), 3) # wonk y-axis
    brush(O - (350, 0), O + (350, 10), 3)  # wonk x-axis

    sine = [(Point(50x, -200sin(x)), Point(50(x + pi/6), -200sin(x + pi/6))) for x in -2pi:pi/10:2pi]

    sethue("pink1")
    brush(Point(-300, -200), Point(300, -200), 0.5)

    sethue("mediumorchid")
    foreach(pr -> brush(pr[1], pr[2], twist=5, highhandle=4, strokes=3), sine)
end 800 600 "images/bezier/graph.png"

a graph

Unfortunately you’ve got the imprecise finish without the reassuring and lovable hand-made quirks. Like those imitation hand-writing fonts, it’s presenting the illusion of manual labour.

@png begin
    fontsize(30)
    fontface("SnellRoundhand-Script")
    text("pint of milk   cat food   stamps   something for dinner? ", halign=:center)
end 800 100 "images/bezier/shoppinglist.png"

a shopping list

A fist full of brushes

But it’s always fun to explore an idea to see where it leads:

@png begin
    t = Table(5, 5, 130, 100)
    l = t.colwidths[1]
    fontsize(4)
    for (pos, n) in t
        weight = rand(0:2:20)
        strokes = rand(1:20)
        minwidth=rand(0.0:0.1:0.3)
        maxwidth=rand(minwidth+0.1:0.1:minwidth+0.3)
        twist = rand(-10:10)
        randomopacity = rand(Bool)
        lowhandle = rand(0.0:0.1:0.5)
        highhandle = rand(lowhandle:0.1:lowhandle+.5)
        tidystart = rand(Bool)
        @layer begin
            translate(pos)
            rotate(-0.2)
            randomhue()
            brush(O - (l/2, 0), O + (l/2, 0), weight,
                strokes=strokes,
                minwidth=minwidth,
                maxwidth=maxwidth,
                lowhandle=lowhandle,
                highhandle=highhandle,
                twist=twist,
                randomopacity=randomopacity,
                tidystart=tidystart)
        end

        textwrap("weight: $weight
            strokes: $strokes
            minwidth: $minwidth
            maxwidth $maxwidth
            twist: $twist
            randomopacity: $randomopacity
            tidystart :$tidystart", t.colwidths[1] - 30, pos + (-l/2, 30))
    end

    sethue("black")
    fontsize(30)
    fontface("Elephant")
    text("brush strokes", boxtop(BoundingBox()) + (0, 40), halign=:center)
end 800 600 "images/bezier/brushes.png"

a catalog of brushes

The line quality can make for simple painterly graphics, good for the occasional Bob Ross painting:

using ColorSchemes
function drawsunset()
    colscheme= ColorSchemes.sunset
    @png begin
        background("burlywood4")
        sethue("orange")
        circle(O + (0, 60), 60, :fill)
        setopacity(0.5)
        for y in -200:20:180
            sethue(get(colscheme, rescale(y, -200, 180, 1, 0)))
            for w in 10:5:20
                brush(Point(-350, y), Point(350, y), w,
                    strokes = 1,
                    twist = 5,
                    minwidth=0.02,
                    maxwidth=0.04,
                    randomopacity=true
                )
            end
        end
        fontsize(30)
        fontface("Pique-Black")
        text("Sunset over Luxor", Point(0, 230), halign=:center)
    end 800 500 "images/bezier/sunset.png"
end
drawsunset()

sunset over luxor

Each time you evaluate this the result is slightly different, yet always the same. Perhaps the next one will be better — wait, no, perhaps I preferred the previous one…

And because we started in France, let’s paint some graffiti:

function drawfrenchflag()
    w, h = 700, 400
    @png begin
        background("azure")
        setopacity(0.85)
        setline(0.5)
        colscheme = [Colors.RGB.(sethue("royalblue1")...),
        Colors.RGB.(sethue("white")...),
        Colors.RGB.(sethue("red")...)]
        for x in -w/2:40:w/2
            # reduce x to 1, 2, or 3, then to [0-1] to get red/white/blue
            col = rescale(div(x+w/2, w/3) + 1, 1, 3, 0, 1)
            sethue(get(colscheme, col))
            for i in 5:10:50
                brush(Point(x, -h/2), Point(x, h/2),
                    i,
                    strokes=2,
                    twist=10,
                    highhandle=1)
            end
        end
        b = box(O, 675, 375, vertices=true)
        for i in 1:length(b)
            sethue("purple")
            brush(b[i], b[mod1(i + 1, length(b))], 30,
                strokes=20,
                minwidth=0.01,
                maxwidth=0.05)
        end
        fontface("MistralStd")
        setopacity(0.7)
        sethue("purple")
        fontsize(100)
        text("Vive la France!", O + (10, 80), halign=:center, angle=-pi/10)
    end 800 500 "images/bezier/frenchflag.png"
end

drawfrenchflag()

drapeau francais

…and then speed off in our curvaceous old Citroen DS and drive to the France/Swiss border, where CERN are busy recreating the big bang:

using ColorSchemes
function unpetitbang()
    @png begin
        background("gray10")

        # rays
        for i in 1:300
            sethue(get(ColorSchemes.inferno, rand()))
            rotate(rand())
            brush(O + (rand(30:500), 0), O + (rand(501:1000), 0), 1,
                strokes=5,
                twist=rand(-1:0.1:1),
                lowhandle  = 0.0,
                highhandle = 0.5
                )
            l = rand(100:600)
            l1 = l + rand(5:50)
            sethue("white")
            brush(O + (l, 0), O + (l1, 0), 1,
                strokes=10,
                twist=rand(-1:0.1:1),
                lowhandle  = 0.0,
                highhandle = 1
                )
        end

        # plasmas
        for i in 1:20
            sethue(get(ColorSchemes.plasma, rand()))
            setopacity(rand(0.25:.05:0.75))
            rotate(rand())
            brush(O - (20, 0), O + (20, 0), 20,
                maxwidth = 10,
                highhandle = 10)
        end

        # orbits
        sethue("ivory")
        setline(0.5)
        for i in 1:10
            rotate(rand())
            brush(O - (0, 0), O + (50, 0), 10,
                maxwidth = 20,
                highhandle = 20,
                action=:stroke)
        end

    end 1000 1000 "images/bezier/unpetitbang.png"
end
unpetitbang()

un petit bang

[2018-06-20]

Footnotes

Don Lancaster

Don Lancaster is the totally awesome dude who was active in the very early days of personal computing, and probably knows more about PostScript than most of the current Adobe Systems employees put together. Travel back in time by visiting his wacky website at http://www.tinaja.com.

Graphic formats

All the images in this post are in PNG, but they looked better in the vector-based SVG format. However, there’s an annoying ‘bug’ in Jupyter/IJulia/IPython involving text in SVG images created by Cairo. What happens is that Cairo tries to be smart and stores text in XML symbols, suitable for re-use. A good idea, but unfortunately they’re stored in the notebook’s ‘global XML scope’, and so later cells accidentally pick up symbol definitions from earlier cells and re-use them, even though that’s not always what you want. A solution would be to somehow encapsulate the SVG image in a cell to prevent the definitions leaking. I don’t know how to do that yet, but there’s an open issue if you can help me find a workround…

This page was generated using Literate.jl.