## 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 practise is commonly referred to as code art, where art is built with code.

The first four posts will aim at defining a function to create a geometric pattern, a spirograph pattern in circle, inspired by this video, with Python libraries.

1. Problem definition and visualization basics in Part1.
2. Line definition through two points and intersection of two lines in Part2.
3. Polygon vertex detection and postprocessing of intersection points Part3.
4. Polygon drawing, colour scheme definition, global function implementation and final drawings in Part4.

We want to make the two previously-developed methods faster by handling the 4-tuple of 4 points that should generate a polygon at once. We need to transform the point arrays into 2 dimensional arrays with as many rows as the number of points to handle simultaneously.

## 2. Vectorized parametric method for polygon vertex detection

We take the points A and B, repeat them 4 times over the vertical axis (row-wise) and assign them to AAs and CCs. I tend to use s at the end of a variable naming to highlight the plurality of the data. Since we want the four lines bordering the polygon, the four rows are listed as:

1. (A, Pka1, B, Pkb1), first segment connects A to first $P_k$, while the second segment connects B to its first $P_k$,
2. (A, Pka1, B, Pkb2), first segment connects A to first $P_k$, while the second segment connects B to its second $P_k$,
3. (A, Pka2, B, Pkb1), first segment connects A to second $P_k$, while the second segment connects B to its first $P_k$,
4. (A, Pka2, B, Pkb2), first segment connects A to second $P_k$, while the second segment connects B to its second $P_k$.

That’s why BBs requires to repeat and DDs to tiles twice vertically.

AAs = np.tile(np.array(Pa), (4,1))
CCs = np.tile(np.array(Pb), (4,1))
print('The single point A {} becomes the 2D array\n{}'.format(Pa, AAs))

The single point A (1.0, 0.0) becomes the 2D array
[[1. 0.]
[1. 0.]
[1. 0.]
[1. 0.]]

idxA, idxB = 12, 18
Pka = (xPs[idxA:idxA+2], yPs[idxA:idxA+2])
BBs = np.repeat(np.vstack(Pka).T, 2, axis=0)
print('The two adjacent points connected to A\n{}\nbecome the 2D array\n{}'.format(Pka, BBs))

The two adjacent points connected to A
(array([-0.80901699, -0.58778525]), array([-0.58778525, -0.80901699]))
become the 2D array
[[-0.80901699 -0.58778525]
[-0.80901699 -0.58778525]
[-0.58778525 -0.80901699]
[-0.58778525 -0.80901699]]

Pkb = (xPs[idxB:idxB+2], yPs[idxB:idxB+2])
DDs = np.tile(np.vstack(Pkb).T, (2,1))
print('The two adjacent points connected to B\n{}\nbecome the 2D array\n{}'.format(Pkb, DDs))

The two adjacent points connected to B
(array([0.80901699, 0.95105652]), array([-0.58778525, -0.30901699]))
become the 2D array
[[ 0.80901699 -0.58778525]
[ 0.95105652 -0.30901699]
[ 0.80901699 -0.58778525]
[ 0.95105652 -0.30901699]]


We build a temporary function to get the points B and D linked to the two generators, given the two indexes.

def index2points(idxA, idxB):
Pka = (xPs[idxA:idxA+2], yPs[idxA:idxA+2])
BBs = np.repeat(np.vstack(Pka).T, 2, axis=0)
Pkb = (xPs[idxB:idxB+2], yPs[idxB:idxB+2])
DDs = np.tile(np.vstack(Pkb).T, (2,1))
return BBs, DDs


Now we need to treat everything as a 2D array where every row potentially gives one of the 4 polygon vertexes. The first vector $\vec{AB}$ is the simple difference between BBs and AAs. Each cross product of two 2D arrays is a 1D array. The cross product of two 2D vectors in theory is a 3D vector, but in practise is treated as a scalar since the first two elements are always 0. Therefore the 1D array coming out of the cross product is the simple list of the four pairs of two 2D vectors (or segments).

ABs = BBs-AAs
cp1 = np.cross(AAs-CCs, DDs-CCs)
cp2 = np.cross(DDs-CCs, ABs)
print('The shape of two 2D arrays cross products is indeed {} and {}'.format(cp1.shape, cp2.shape))

The shape of two 2D arrays cross products is indeed (4,) and (4,)


We need to make sure cp2 does not contain 0 values (that happens when $\vec{CD}$ and $\vec{AB}$ are parallel) to prevent either NaN or $\infty$ values as a result of the division operation. So first we keep note of such events and then we simply replace every 0 with whatever else (say, 1).

paralSeg = cp2==0
cp2[paralSeg] = 1


To apply element-wise product between the 2D array ABs and the 1D array got from the ratio of two cross products, by exploiting Numpy broadcasting, we need to reshape the latter array to a 2D with .reshape(-1, 1). Finally the (4, 2) array is split into two (4,) arrays, which represent the $(x, y)$ coordinates of the four intersecting points.

print('Shape of ratio of two cross products is {}'.format((cp1/cp2).shape))

Shape of ratio of two cross products is (4,)

PIs = ABs*(cp1/cp2).reshape(-1,1)+AAs
xIs, yIs = np.split(PIs, 2, axis=1)


We show here the four lines defining one of the possible polygons and the four intersecting points (green dots). We draw the points only if the lines are coincident (paralSeg==0).

figSize = 10
fourLines(idxA, idxB)
plt.plot(xIs[paralSeg==0], yIs[paralSeg==0], ls='', marker='o', markersize=10, color="green");


We check the code for a case of parallel/collinear segments.

idxA, idxB = 6, 0
BBs, DDs = index2points(idxA, idxB)

ABs = BBs-AAs
cp1 = np.cross(AAs-CCs, DDs-CCs)
cp2 = np.cross(DDs-CCs, ABs)
are2segParallel = np.abs(cp2)<1e-5
cp2[are2segParallel] = 1

PIs = ABs*(cp1/cp2).reshape(-1,1)+AAs
xIs, yIs = np.split(PIs, 2, axis=1)

figSize = 10
fourLines(idxA, idxB)
plt.plot(xIs[paralSeg==0], yIs[paralSeg==0], ls='', marker='o', markersize=10, color="green");


## 3. Vectorized homogeneous method for polygon vertex detection

We have got the 4 tuples of 4 points A, B, C and D stored with same structure in arrays AAs, BBs, CCs and DDs.

We transform each of the four points into homogeneous coordinates, namely $(x_a, y_a)$ in AAs becomes $(x_a, y_a, 1)$ in AHs, and so on.

idxA, idxB = 12, 18
BBs, DDs = index2points(idxA, idxB)

AHs = np.hstack((AAs, np.ones((4, 1))))
BHs = np.hstack((BBs, np.ones((4, 1))))
CHs = np.hstack((CCs, np.ones((4, 1))))
DHs = np.hstack((DDs, np.ones((4, 1))))
print('Each of these 2D arrays has shape equal to {}'.format(AHs.shape))

Each of these 2D arrays has shape equal to (4, 3)


L1 is the coordinates of the line through points A and B, as a cross product of the homogeneous coordinates of those points. The shape does not change.

L1 = np.cross(AHs, BHs)
L2 = np.cross(CHs, DHs)


PIHs is the homogeneous coordinates for intersecting point of the two lines. Every row is the intersection between one of the four-line pair combinations. We horizontally (axis=1) split the (4, 3) arrays into the three homogeneous coordinate elements, xIHs, yIHs and zIHs, each being a 1D array with shape (4,).

PIHs = np.cross(L1, L2) # homogeneous coordinates for intersecting point
xIHs, yIHs, zIHs = np.split(PIHs, 3, axis=1)


We check that the 3D scaling factor zIHs is not 0, otherwise the lines are parallel or collinear. We finally convert homogeneous into Euclidean coordinates by dividing by zIHs.

paralSeg = np.abs(zIHs)<1e-5
zIHs[paralSeg] = 1
xIs, yIs = xIHs/zIHs, yIHs/zIHs


We draw the points only if the lines are coincident (paralSeg==0).

figSize = 10
fourLines(idxA, idxB)
plt.plot(xIs[paralSeg==0], yIs[paralSeg==0], ls='', marker='o', markersize=10, color="green");


We check the code for a case of parallel/collinear segments.

idxA, idxB = 6, 0
BBs, DDs = index2points(idxA, idxB)

BHs = np.hstack((BBs, np.ones((4, 1))))
DHs = np.hstack((DDs, np.ones((4, 1))))
L1 = np.cross(AHs, BHs)
L2 = np.cross(CHs, DHs)
PIHs = np.cross(L1, L2) # homogeneous coordinates for intersecting point
xIHs, yIHs, zIHs = np.split(PIHs, 3, axis=1)
paralSeg = np.abs(zIHs)<1e-5
zIHs[paralSeg] = 1
xIs, yIs = xIHs/zIHs, yIHs/zIHs

figSize = 10
fourLines(idxA, idxB)
plt.plot(xIs[paralSeg==0], yIs[paralSeg==0], ls='', marker='o', markersize=10, color="green");


## 4. Postprocess the set of intersecting points

We need to process the list of coordinates for the detected intersecting points.

Here the main steps to perform:

1. We first need to drop any duplicates, which means points that coincide with each other. Despite Matplotlib can handle polygons with duplicated coordinates, it is a good practise to remove those points.
2. Whatever point lying outside the circle has to be removed.
3. The array of point coordinates needs to be sorted in a (anti-)clockwise fashion to make sure the patch Matplotlib has to fill with colour is convex. A simple visual example will come soon.

### 4.1 Remove duplicate points

Let’s take a case where we end up with duplicate intersecting points (idxA=18 and idx=19).

idxA, idxB = 18, 19
BBs, DDs = index2points(idxA, idxB)

BHs = np.hstack((BBs, np.ones((4, 1))))
DHs = np.hstack((DDs, np.ones((4, 1))))
L1 = np.cross(AHs, BHs)
L2 = np.cross(CHs, DHs)
PIHs = np.cross(L1, L2) # homogeneous coordinates for intersecting point
xIHs, yIHs, zIHs = np.split(PIHs, 3, axis=1)
paralSeg = np.abs(zIHs)<1e-5
zIHs[paralSeg] = 1
xIs, yIs = xIHs/zIHs, yIHs/zIHs

figSize = 8
fourLines(idxA, idxB)
plt.plot(xIs[paralSeg==0], yIs[paralSeg==0], ls='', marker='o', markersize=10, color="green");


You can realize the point $(1, 0)$ is repeated twice, placed at second and last rows. However we cannot simply feed that to the Numpy unique function since the two points are not numerically identical. First those points need to be rounded up to the 6-th decimal digits with the around method and then fed to the unique function.

polygonalPoints = np.hstack((xIs, yIs))
print('Full set of coordinates identified for the 4 pairs of intersecting lines:\n{}'.format(polygonalPoints))

Full set of coordinates identified for the 4 pairs of intersecting lines:
[[ 9.09422702e-01 -2.78768258e-01]
[ 1.00000000e+00 -1.53094054e-16]
[ 9.51056516e-01 -3.09016994e-01]
[ 1.00000000e+00 -2.09377778e-16]]

polygonalPoints1 = np.unique(np.around(polygonalPoints, decimals=6), axis=0)
print('Set of unique rounded coordinates identified for the 4 pairs of intersecting lines:\n{}'.format(polygonalPoints1))

Set of unique rounded coordinates identified for the 4 pairs of intersecting lines:
[[ 0.909423 -0.278768]
[ 0.951057 -0.309017]
[ 1.       -0.      ]]


### 4.2 Remove points outside the circle

In the original drawing every polygon lies within the circle only, so we need to remove every point of the intersecting set exceeding it. The logic is quite straightforward. First the squared distance of the point from the circle centre $(0, 0)$ is obtained by squaring the coordinates element-wise, summing along the columns and finally taking the square root to get a linear distance. We round the distance to the third digit and check it is less than the circle radius radius. The intersecting point set is reduced to such points that satisfy this check only. In Numpy the syntax to select some rows of array AA that meet some criteria cc is AA[cc==True, :].

distance = np.sqrt(np.sum(polygonalPoints1**2, axis=1))
polygonalPoints2 = polygonalPoints1[criteria, :]


### 4.3 Sort the set of intersecting points anticlockwise

We generate 4 2D points as vertexes of a polygon and draw it in purple on the LHS, making sure that it is concave (non-convex).

A swapped version of the same set (2 and 3 rows) is then shifted by 5 units rightward and depicted in yellow. Clearly the two shapes are quite different, according to the specific order of the points along the array rows. We need to sort the set in some way. We take the anticlockwise sense.

polygCoord = np.array([[2,-1], [0,1], [1,1.5], [-1,-1.5]])
polygCoord1 = polygCoord+np.array([5, 0])
swapIdx = [1, 2]
polygCoord1[swapIdx] = polygCoord1[swapIdx[::-1]]

plt.figure(figsize=(10, 5))
patches = [Polygon(polygCoord, True), Polygon(polygCoord1, True)]
pc = PatchCollection(patches, alpha=.3)
pc.set_array(np.array([0, 1]))

ax = plt.gca()

ax.get_xaxis().set_ticks([])
ax.get_yaxis().set_ticks([])
ax.axis('off')
ax.axis('equal');


We now complete the logic to sort anticlockwise.

If there are still 4 intersecting points after applying the two previous steps, we need to sort those points anticlockwise.

Here it comes the polar coordinate system at hand one more time.

We first calculate the centroid C of the polygon as the geometric mean of the 4 points. It is the row-wise (axis=0) mean of point coordinates set (red point the below chart).

Then the angle $\theta_j$ is formed between the segment connecting the $P_j$ and the centroid C and the x axis. That angle gives the criterion to sort the 4 points. Since angle increases anticlockwise by convention in polar system, points are sorted accordingly, but it has no effect on the polygon visualization in Matplotlib.

The angle between one segment $\vec{CP_j}$ and the horizontal line is given by the arctangent function as:

$$\theta_j = \arctan{\frac{y_{P_j}-y_C}{x_{P_j}-x_C}}$$

This function is implemented both in Numpy and in the math library as atan2(y, x).

Since we need to sort the 2D array of point coordinates $(x,y)$ with respect to a third temporary variable, $\theta$, the Python sorting functionality is required due to a lack of similar package in Numpy. Indeed, one can here use the key attribute to the sort function to give the sorting criterion. That’s why we use the atan2 from the math library.

The new polygon is shifted by 5 units leftward and depicted in green.

centroid = np.mean(polygCoord, axis=0)
polygCoord_ = list(polygCoord)
polygCoord_.sort(key=lambda pnt: math.atan2(pnt[1]-centroid[1], pnt[0]-centroid[0]))
polygCoord2 = np.array(polygCoord_)-np.array([5, 0])
patches.insert(1, Polygon(polygCoord2, True))

plt.figure(figsize=(20, 5))
plt.plot(centroid[0], centroid[1], ls='', marker='o', markersize=10, color="red");
pc = PatchCollection(patches, alpha=.3)
pc.set_array(np.array([0, 1, 2]))
ax = plt.gca()