1. Introduction

This post belongs to a new series of posts where I intend to face the challenge of drawing both static and dynamic pictures with a programming language. This practice is commonly referred to as code art, where art is built with code.

This post and the next one will aim at defining a function to create a geometric pattern, a spirograph pattern in a rectangle, inspired by this video, with Python libraries.

We will follow these steps:

  1. Grid definition, parametric equation of the rectangle in Part1.
  2. Polygon vertex detection, postprocessing of intersection points, colour scheme definition, global function implementation and final drawings in Part2.

2. Description

We develop the general function which is responsible to draw the spirograph pattern art in a rectangle.

We define the same pattern as we did in the previous post-series (from Part1 to Part4). The only difference is that the pattern is based on a rectangular grid, rather than a circle.

The hyperparameters are:

  1. Nh, the number of points over the rectangle height.
  2. height, the rectangle height.
  3. width, the rectangle width.
  4. genLocs, the relative distance of the two generators, A and B, from the leftmost rectangle side. Both generators lie on the horizontal sides. The A generator lies on the bottom side, while B lies on the top side.
  5. kind, the kind of color scheme, which can be one among generator-distance based genDist, fix-color fix-(green|red|blue|gray|black), polygon-area-based area and origin-distance-based origiDist schemes.
  6. figSize, the figure size
  7. palette, the colour palette from Matplotlib.
  8. patchAlpha, the polygon colour transparency.

3. Rectangular grid

Let a segment be a fraction of the rectangle height ($h$), whose length Lseg ($L$) is defined as follows:

$$ L = \frac{h}{N_h} $$

where $N_h$ is the number of points Nh.

For the sake of generality, we let the user define the width as an independent parameter. We therefore need to estimate the actual width ($w$) once again, as the product of the segment length and the least integer greater than or equal to the ratio between the original width and the segment length itself:

$$ w = L\cdot ceil\Big(\frac{w}{L}\Big) $$

Recall the ceiling function definition, as the least integer greater than or equal to the input.

We also ask the user to specify where to place the two generators, A and B. The parameter genLocs ($\delta$) is a tuple of two values ranging from 0 to 1. The coordinates of the two generators are thus defined:

$$ A = (x, y) = (w\cdot \delta_0, 0) $$

$$ B = (x, y) = (w\cdot \delta_1, h) $$

where the expression $\delta_j$ is the j-th element of the tuple.

However, the two generators have to be placed exactly on any of the grid nodes of the two horizontal sides. We introduce the function coord2index to convert the general coordinates of the two generators into grid coordinates. The nearest node is selected with the round function as the new generator location.

coord2index = lambda xx, Lseg: np.round(xx/Lseg)

Figure 1 explains how we can discretize the rectangle with a grid. In this case, the height is split into 4 segments, or units, while the width is the 8-unit long.


Figure 1 - Rectangle grid and indexing scheme

We use the following snippet of code to determine the actual width, the indexes and the coordinates of the two generators.

width, height, Nh = 8, 4, 4
genLocs=(3/8, 5/8)

Lseg = height/Nh
Nw = coord2index(width, Lseg)
Ntot = int(2*(Nw + Nh))
width = Nw*Lseg

kA = coord2index(genLocs[0]*width, Lseg)
kB_ = coord2index(genLocs[1]*width, Lseg)
kB = 2*Nw + Nh - kB_
AA = np.r_[kA*Lseg, 0]
BB = np.r_[kB_*Lseg, height]
AAs = np.tile(AA, (4,1))
BBs = np.tile(BB, (4,1))
idxAB = kB - kA - 1
print('The coordinates of the two generators A and B are {} and {}, respectively'.format(AA, BB))
print('The index distance between these two generators is {}'.format(idxAB))
The coordinates of the two generators A and B are [3. 0.] and [5. 4.], respectively
The index distance between these two generators is 11.0

According to the distance of the two generators from the left side of the rectangle, we end up with any of these two different cases (and shapes):

  1. A is on the left- or right-hand side of the vertical line crossing B, so the central line is inclined.
  2. They are equally distant from the leftmost vertical side of the rectangle, so the central line is vertical.

4. Converting grid index into time

Let’s start the actual process of generating polygons by connecting A and B to any two consecutive nodes on the grid. This means that we have two pairs of lines that are going to potentially intersect and give birth to several polygons.

The process is managed by counter-clockwise indexing the nodes on the rectangular grid from the bottom-leftmost node, $(0,0)$. We want to have two increasing indexes, dka and dkb, to progressively pick a pair of consecutive nodes to connect to the generator, A and B, respectively, within two nested for-loops, one per generator, as you can see in the main function code drawingSpirograph.

We need to find a way to convert the monotonic increasing index into a specific position within the four sides of the rectangle and then finally got the coordinates. When we have the coordinates of the four lines, we can execute the function that returns the coordinates of the intersecting points for the four lines.

With respect to Figure 1, we observe that k goes from 0, which is the origin, to $N_w=8$, which is in general the bottom rightmost corner, to $N_w+N_h=12$, which is the top rightmost point, until it comes back to the origin when it equals $2\cdot(N_w+N_h)=24$. $N_w$ is the number of intervals of units for the width, while $N_h$ is the number of intervals of units for the height.

The index specifies where the point to select is and in which side of the rectangle it is. But we need a second approach to identify what’s the correlation between the position and the actual coordinates.

To visualize the following approach, we just imagine to draw the rectangle with a pencil, one side at a time. Each side is drawn in one unit of time. We start from the origin, we draw the bottom line in 1 second to end up being in point $(w,0)$. After one second more, we complete the second side to reach $(w, h)$. We draw the top line during the third second and close the shape in the fourth unit of time.

That means that, since the time to draw every side is constant but its length is in principle different, we have a different speed of drawing. We are therefore faster when we draw the horizontal lines and slower for vertical ones, which is pretty acceptable because, to some extent, we only care about the final shape but we don’t consider any kinematics/dynamics aspects in this particular occasion.

The below chart correlates the grid index to the time at which the indexed point is drawn. The below function implements this conversion from k to time t. It is a piecewise-linear function combined with a powerful interpolation function in numpy. We can just give the input k, get rid of extra-lap gap Ntot with the modulus operator % and interpolate the function to get the time.


Figure 2 - Rectangle index-to-time conversion

def rectIndex2Time(kk, Nw, Nh, Ntot):
    kk = kk%(Ntot)
    return np.interp(kk, [0, Nw, Nw+Nh, 2*Nw+Nh, Ntot], list(range(5)))

5. Converting time into coordinates

We need a function that specifies how the position of the point with respect to the x-axis changes over time from 0 to 4 seconds. This would give the parametric equation of the first coordinate of the point. Similarly, we introduce an additional function for the y-axis. Altogether, we have a system of parametric equations, that is, the parametric equation of the rectangle.

Let’s just summarize the idea! From the grid, we have a specific position. We are able to convert the grid index to the corresponding time at which we pass by while drawing the rectangle. We need a second step to convert this time variable into a cartesian coordinate variable.

While we imagine drawing the rectangle, we observe how the $x$ coordinate changes. It first increases from 0 to $w$ over the first second, it is constant for one second, it decreases down to zero in the third second and it finally doesn’t change any more.

The $y$ coordinate is 0 for one second, increases to $h$ in one second, it maintains this value for one second and it drops to 0 during the last second.

How can we define a function that converts time to the two coordinates, x and y?

We have two different approaches. The former is what we use for the index-to-time conversion, i.e., a piecewise-linear function combined with a powerful interpolation function in numpy. This time-to-coordinate conversion requires two piecewise-linear functions, one per coordinate, x and y. But just to get a broad scope, we can use an alternative method that combines the equation of every segment and then uses a combination of the min and max operators, to end up having the overall piecewise-linear function.


Figure 3 - Rectangle time-to-coordinate conversion, x-axis only

With reference to Figure 3, let’s start realizing that the first segment equation is the line $x=t$ (green), the second segment is $x=1$ (yellow), the third one is $3-t$ (blue) and the last one is $x=0$.

How to combine these four equations into a single one? If we take the minimum between the green and blue lines we have a kind of roof, $x = min(t, 3-t)$. The actual roof is capped, it has a trapezoid shape. So we take again the minimum between the roof and the yellow line, $x = min(1, min(t, 3-t))$. The x coordinate would however be negative for $t>3$, so we use the max operator as:

$$ x(t) = max\big(0, min(1, min(t, 3-t)\big) $$


Figure 4 - Rectangle time-to-coordinate conversion, x- and y-axis

To determine the y coordinate parametric equation, we can simply apply a trick (see Figure 4). Since its shape is exactly identical to x but shifted by one unit of time, we can take the same function $f(t)$ that describes x and gives the shifted time $t-1$ as input.

The overall system is as compact as:

$$ x(t) = f(t) $$

$$ y(t) = f(t-1) $$

The below function implements this conversion from time to coordinates, xx and yy.

def rectTime2Coord(time, width, height):
    roof = lambda time: np.minimum(time, 3-time)
    xx = width*np.maximum(np.zeros_like(time), np.minimum(roof(time), np.ones_like(time)))
    yy = height*np.maximum(np.zeros_like(time), np.minimum(roof(time-1), np.ones_like(time)))
    return np.array((xx, yy)).T

6. Full conversion from index to coordinates

From the indexes ka and kb, which are the indexes of the two generators, we get the coordinates of the two pairs of adjacent points that are then connected to the two generators to design the polygon.

We first get the temporal instant for the C-points pair, whose indexes are dka+1 and dka+2 away from A, by using rectIndex2Time. The second step is to get the coordinates of the C-points pair with rectTime2Coord. We repeat these two steps for the D-points pair. To get the full combination of the two pairs of C and D, we repeat the C coordinates and we tile the D coordinates. This logic is in rectIndex2coords.

dka, dkb = 3, 5
tCs = rectIndex2Time(kA + (np.r_[dka, dka+1]+1), Nw, Nh, Ntot)
Pcs = rectTime2Coord(tCs, width, height)
print('The temporal instants for the C-points pair are {} and {}'.format(*tCs))
print('The coordinates are {} and {}'.format(*Pcs))
The temporal instants for the C-points pair are 0.875 and 1.0
The coordinates are [7. 0.] and [8. 0.]
def rectIndex2coords(kA, kB, width, height, dka, dkb, Nw, Nh, Ntot):
    tCs = rectIndex2Time(kA + (np.r_[dka, dka+1]+1), Nw, Nh, Ntot)
    Pcs = rectTime2Coord(tCs, width, height)
    CCs = np.repeat(Pcs, 2, axis=0)
    tDs = rectIndex2Time(kB + (np.r_[dkb, dkb+1]+1), Nw, Nh, Ntot)
    Pds = rectTime2Coord(tDs, width, height)
    DDs = np.tile(Pds, (2,1))
    return CCs, DDs
CCs, DDs = rectIndex2coords(kA, kB, width, height, dka, dkb, Nw, Nh, Ntot)

The set of the four points’ coordinates that gives the vertexes of the polygon is printed as a Pandas dataframe, where column names combine the double point name (from A to D) with one of the two axes.

import pandas as pd
colNames = [point*2+'-'+axisType for point in 'ABCD' for axisType in 'xy']
pd.DataFrame(np.c_[AAs, BBs, CCs, DDs], columns=colNames)
AA-x AA-y BB-x BB-y CC-x CC-y DD-x DD-y
0 3.0 0.0 5.0 4.0 7.0 0.0 0.0 3.0
1 3.0 0.0 5.0 4.0 7.0 0.0 0.0 2.0
2 3.0 0.0 5.0 4.0 8.0 0.0 0.0 3.0
3 3.0 0.0 5.0 4.0 8.0 0.0 0.0 2.0