A surprisingly wide variety of interesting-looking vector line graphics can be generated with a fairly simple technique.

First, we take a bunch of lines. These lines work best if they represent a pattern with some regularity, so I call them the **texture**. Then, we pick a **surface**, which is essentially an elevation map represented by a matrix of heights. It can be generated mathematically or taken from real data. Next, we **map** the texture onto the surface, and capture the heights along all each line to add a third dimension to those lines. Finally, **project** the result to convert the 3D lines back into 2D lines. If we're lucky, these lines will be more interesting to look at than the texture we started with.

We'll use penkit's `show_layer`

function for displaying plots within the notebook, so if the following code throws an `ImportError`

, try running `pip install penkit`

in a terminal.

In [1]:

```
%matplotlib inline
import matplotlib.pyplot as plt
from penkit.preview import show_layer
import numpy as np
```

Let's begin by generating a texture. As I mentioned, a texture is just a set of lines.

For the first texture we'll take inspiration from the band Joy Division's famous Unknown Pleasures album cover. Joy Division-inspired plots have become a genre in their own right, sometimes called joy plots

The texture in this case is simply a bunch of evenly-spaced horizontal lines. To simplify the math later on, it is helpful to constrain texture coordinates such that both dimensions are in the range `[0, 1]`

. If you try generating your own (and I encourage you to!), make sure they are within this area (you might find `penkit.textures.util.fit_texture`

helpful).

In [2]:

```
def make_joy_texture(num_lines=10, resolution=50):
# np.meshgrid is a handy way to generate a grid of points. It
# returns a pair of matrices, which we will flatten into arrays.
# For the x-coordinates, we put a nan value at the end so that when
# we flatten them there is a separater between each horizontal line.
x, y = np.meshgrid(
np.hstack([np.linspace(0, 1, resolution), np.nan]),
np.linspace(0, 1, num_lines),
)
# For coordinates where the x value is nan, set the y value to nan
# as well. nan coordinates represent breaks in the path, indicating
# here that the pen should be raised between each horizontal line.
y[np.isnan(x)] = np.nan
return x.flatten(), y.flatten()
```

`nan`

("not a number") indicator, which represents a gap in the path. All coordinates are absolute, i.e. vectors from the origin rather than relative to the previous point. If you would like a more gradual introduction to this way of representing lines, see my Fractal Generation with L-Systems tutorial.

In [3]:

```
joy_texture = make_joy_texture(6)
```

We can use `show_layer`

to preview the texture:

In [4]:

```
show_layer(joy_texture)
```

Out[4]:

`num_lines`

. You might have also noticed the optional `resolution`

parameter and wondered why it matters. The lines above are actually made up of a bunch of line segments, and the number of these segments is the **resolution**. We can use `matplotlib`

to make it clear where the line segments are.

In [5]:

```
plt.plot(*make_joy_texture(10, 5), marker='.');
```

`make_joy_texture`

to check your understanding. The resolution doesn't make a difference when we plot the texture, but soon we will be bending the texture into curves. The higher the resolution, the more accurate the curves will be.

In [6]:

```
from scipy.ndimage.filters import gaussian_filter
# This step isn't necessary, but it ensures that the same "random" plot
# will be generated every time the tutorial is run.
np.random.seed(int('banana', 36))
NOISE_BLUR = 40
noise_surface = gaussian_filter(np.random.normal(size=(500, 500)), NOISE_BLUR)
```

`matplotlib`

, but we need to use `imshow`

rather than `plot`

. Lighter pixels represent points with higher elevations.

In [7]:

```
plt.imshow(noise_surface, cmap='gray');
```

Another way to generate a surface is to run a mathematical function over the x and y coordinates. To do this easily in `numpy`

, we can generate a pair of matrices:

In [8]:

```
SURFACE_SIZE = 500
grid = np.meshgrid(
np.linspace(0, 1, SURFACE_SIZE),
np.linspace(0, 1, SURFACE_SIZE)
)
```

In [9]:

```
_, (p1, p2) = plt.subplots(1, 2)
p1.imshow(grid[0], cmap='gray')
p2.imshow(grid[1], cmap='gray');
```

Standard graphing functions, like `sin`

, can be generalized into 3D by converting the (x, y) values into a scalar and calling `sin()`

on *that*. The natural transformation into a scalar is the *norm*, i.e. the distance from the origin.

Note that the x and y values are scaled between 0 and 1, rather than being the raw pixel coordinates.

In [10]:

```
hole_surface = np.sin(
np.linalg.norm((np.array(grid) - 0.5) * np.pi, axis=0))
```

In [11]:

```
plt.imshow(hole_surface, cmap='gray');
```

Another fun math function involves taking the product of the sine of both axis. It comes out looking like a grid of bubbles, so I call it the `bubble_surface`

.

In [12]:

```
NUM_BUBBLES_PER_SIDE = 3
bubble_surface = (
np.sin((grid[0] - 0.5) * NUM_BUBBLES_PER_SIDE * np.pi) *
np.sin((grid[1] - 0.5) * NUM_BUBBLES_PER_SIDE * np.pi))
```

In [13]:

```
plt.imshow(bubble_surface, cmap='gray');
```

Now that we have a surface, we can **map** it to the texture. First we strech the texture over the surface dimensions. Then, for each point in the texture, we find its associated pixel on the surface and look up the "elevation" at that point. That elevation becomes the Z-value for the texture. We can think of these Z-values, combined with the original X and Y values, as making up a 3D version of the original texture.

In [14]:

```
def texture_map(texture, surface):
texture_x, texture_y = texture
surface_w, surface_h = surface.shape
# First, we convert the points along the texture into integers within
# the bounds of the surface's index. The clipping here will also convert
# the nan values to 0.
index_x = np.clip(np.int32(surface_w * texture_x), 0, surface_w - 1)
index_y = np.clip(np.int32(surface_h * texture_y), 0, surface_h - 1)
# Grab z-values along the texture path. Note that this will include values
# for every point, even if it is nan or had to be clipped to within the
# bounds of the surface, so we have to fix that next.
surface_z = surface[index_x, index_y]
# Set z-values that are either out of bounds or nan in the texture to nan
# in the output.
# Numpy wants to warn about the fact that there are nan values in the
# textures, so we silence warnings about this.
with np.errstate(invalid='ignore'):
surface_z[(texture_x < 0) | (texture_x >= 1) |
(texture_y < 0) | (texture_y >= 1)] = np.nan
return surface_z
```

**project** it back to 2D. There are different ways to do this, but one way is to preserve the X coordinates, and use a linear blend of the Y and Z coordinates as the Y coordinates. Here's a function that handles both the **mapping** (by calling `texture_map`

) and **projection**.

In [15]:

```
def texture_plot(texture, surface, angle=45, **kwargs):
# Extract the Xs and Ys from the texture
surface_x, surface_y = texture
# Map the texture to get the Zs
surface_z = texture_map(texture, surface.T)
# The projection is as simple as linearly blending the Z and Y
# dimensions. The multiples are calculated from the given
# angle with standard trig.
z_coef = np.sin(np.radians(angle))
y_coef = np.cos(np.radians(angle))
plot = (surface_x, -surface_y * y_coef + surface_z * z_coef)
return show_layer(plot, **kwargs)
```

In [16]:

```
texture_plot(make_joy_texture(20), bubble_surface * 0.1, 40)
```

Out[16]:

In [17]:

```
texture_plot(make_joy_texture(40), noise_surface * 12, 40)
```

Out[17]:

You may have noticed that in the plot above, some of the lines cross each other. This gives the impression that the surface is transparent, as if we can "see through" it to lines behind the surface. In order to make images that look 3D, we can make the surface appear opaque by not drawing lines that would be obscured by the surface.

One way to do this is by "masking" the visible parts of the surface. First, we apply the projection to the surface. Then, we walk over the surface keeping track of the maximum height as we go. Only the parts of the surface where the height of the sloped surface is equal to the (cumulative) maximum are kept.

To understand this, consider a two-dimensional version below. The lines of sight are parallel, because we are using an axonometric projection. In this example, the lighter-colored parts of the surface would be masked out because they are hidden to the observer. If you trace the line left to right, you can see that points along the line are only visible if no point left of them is higher. In other words, they are equal to the cumulative maximum point.

`matplotlib`

has a function `np.maximum.accumulate`

that calculates a cumulative maximum over every slice at once. By comparing the cumulative maximums to the actual elevation, we can determine the visibility of each pixel on the surface. The result of this function is a mask that tells us which pixels of the surface are visible.

In [18]:

```
def make_visible_mask(surface, angle):
s = surface.shape[0]
# Just as in projection, we calculate the Y and Z
# coefficients with sin/cos.
y_coef = np.cos(np.radians(angle))
z_coef = np.sin(np.radians(angle))
# Rotate the surface so that the visibilty mask represents
# the visibility at the desired angle.
projected_surface = (
z_coef * surface -
y_coef * np.expand_dims(np.linspace(0., 1., s), axis=1)
)
# Calculate the cumulative maximum along each cross-section of
# the projected surface. We flip on the input and output because
# we want to accumulate from the bottom of the surface up, rather
# than the top-down. This is because we interpret the bottom of
# the surface as closer to the observer.
projected_surface_max = np.flipud(np.maximum.accumulate(np.flipud(projected_surface)))
# Compare each point on the surface to the cumulative maximum
# along its cross-section.
return projected_surface == projected_surface_max
```

To apply the mask, we can set all values of the surface that are *not* in the mask to `nan`

.

In [19]:

```
def remove_hidden_parts(surface, angle):
surface = np.copy(surface)
surface[~make_visible_mask(surface, angle)] = np.nan
return surface
```

The input and output of this process, as well as the mask used, are shown below.

In [20]:

```
_, (i, m, r) = plt.subplots(1, 3)
# Plot input surface
i.imshow(bubble_surface, cmap='gray')
i.set_title('input')
i.axis('off')
# Plot mask
m.imshow(make_visible_mask(bubble_surface, 10), cmap='gray')
m.set_title('mask')
m.axis('off')
# Plot result
r.imshow(remove_hidden_parts(bubble_surface, 10), cmap='gray')
r.set_title('result')
r.axis('off');
```

Putting it all together, we can apply `remove_hidden_parts`

to a surface and get the result we want.

In [21]:

```
angle=70
texture_plot(make_joy_texture(50, 100), remove_hidden_parts(bubble_surface * 0.1, angle), angle)
```

Out[21]:

In [22]:

```
angle=65
texture_plot(make_joy_texture(80, 100), remove_hidden_parts(
hole_surface * 1.2 + noise_surface + bubble_surface * 0.1, angle), angle)
```

Out[22]: