A Python 3 script that generates art using matplotlib

Question

I am always super interested in both science and art, science and art are two different ways to describe the objective reality, they are two sides of the same coin, in many areas they overlap, and mathematics is the most prominent one of such areas.

So I used the last three days to write a bunch of functions that generate artistic geometrical shapes using mathematics and visualize the shapes using matplotlib.

Previously I had absolutely zero experience with matplotlib, in the past I have used GeoGebra 5 to draw shapes, but doing it manually is really slow and prone to error, and I don't use brushes because I want my images to be precise, so I have decided to draw pictures with code and here I am.

In this gigantic script I have combined:

Hailstone Sequence

Fibonacci Sequence

Prime Factorization

Dice Average Formula

The following equation is used to calculate the average value of a set of dice.

AV = ( (M + 1 )/ 2 ) * N

Where AV is the average dice value
M is the max value of all dice
N is the total number of die

source

Full Saturation Spectrum

Explanations by me: the RGB color model uses one byte for each channel to represent colors, one byte is eight bits, 2 to the power of 8 is 256, so each channel can have a value from 0 to 255, and the total number of colors possible in this model is 16777216 (2 to the power of 24).

I don't know about much about color theory, but in my observation, fully saturated colors have only one or two bytes that are set, and at least one of the bytes is fully set. So the total number of fully saturated colors is 6 * 255 = 1530.

And of course, programming, I have written over 400 lines of code, spent more than 12 hours on this project, expressed math in code, just to create, something with absolutely no practical value whatsoever—ART.

---

Code

from functools import lru_cache
from matplotlib.patches import Arc, Circle, Rectangle
from itertools import combinations
import matplotlib.pyplot as plt
import random
import math
import re

FIBONACCI = [2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233]
DIRECTIONS = [('right', 'up'), ('right', 'down'), ('left', 'down'), ('left', 'up')]
ANGLES = [(90, 180), (0, 90), (-90, 0), (180, 270)]
CCW_ANGLES = [(0, 90), (90, 180), (180, 270), (270, 360)]
cycles = [
    (4, 2, 1),
    (-2, -1),
    (-14, -7, -20, -10, -5),
    (
        -50, -25, -74,
        -37, -110, -55,
        -164, -82, -41,
        -122, -61, -182,
        -91, -272, -136,
        -68, -34, -17
    )
]
loops = {term: len(cycle) for cycle in cycles for term in cycle}
stops = {term: cycle[-1] for cycle in cycles for term in cycle}

@lru_cache(maxsize=None)
def fibonacci(n):
    s = []
    a, b = 0, 1
    for i in range(n):
        s.append(a)
        a, b = b, a + b
    return s

@lru_cache(maxsize=None)
def spectrum_position(n):
    if not isinstance(n, int):
        raise TypeError('`n` should be an integer')
    if n < 0:
        raise ValueError('`n` must be non-negative')
    n %= 1530
    if 0 <= n < 255:
        return f'#ff{n:02x}00'
    elif 255 <= n < 510:
        return f'#{510-n:02x}ff00'
    elif 510 <= n < 765:
        return f'#00ff{n-510:02x}'
    elif 765 <= n < 1020:
        return f'#00{1020-n:02x}ff'
    elif 1020 <= n < 1275:
        return f'#{n-1020:02x}00ff'
    elif 1275 <= n < 1530:
        return f'#ff00{1530-n:02x}'

@lru_cache(maxsize=None)
def factorize(n):
    assert isinstance(n, int)
    assert n >= 2
    factors = []
    while not n % 2:
        factors.append(2)
        n /= 2
    
    for i in range(3, int(n**.5+1), 2):
        while not n % i:
            factors.append(i)
            n /= i
    
    if n > 1:
        factors.append(int(n))
    return factors

@lru_cache(maxsize=None)
def dice_average(n):
    factors = factorize(n)
    if len(factors) == 1:
        return (1, n * 2 - 1)
    else:
        maxf = max(factors)
        factors.remove(maxf)
        return (math.prod(factors), maxf * 2 - 1)

@lru_cache(maxsize=None)
def hailstone(n):
    if not isinstance(n, int):
        raise TypeError('n should be an integer')
    
    if n == 0:
        raise ValueError('n should not be 0')
    
    sequence = [n]
    loop_iter = 0
    looping = False
    while True:
        if not looping and n in loops:
            loop_terms = loops[n]
            looping = True
        
        if looping:
            loop_iter += 1
            if loop_iter == loop_terms:
                break
        
        if n % 2:
            n = 3 * n + 1
        else:
            n /= 2
        
        sequence.append(int(n))
    
    return sequence

def sin(x): return math.sin(math.radians(x))
def cos(x): return math.cos(math.radians(x))
def tan(d): return math.tan(math.radians(d))

@lru_cache(maxsize=None)
def negate(ls): return [-e for e in ls]

@lru_cache(maxsize=None)
def rotate(point, deg):
    x, y = point
    x1 = x * cos(deg) - y * sin(deg)
    y1 = y * cos(deg) + x * sin(deg)
    return (x1, y1)

@lru_cache(maxsize=None)
def square_vertices(init_pos, side, direction):
    x0, y0 = init_pos
    hd, vd = direction
    x1 = x0 + side if hd == 'right' else x0 - side
    y1 = y0 + side if vd == 'up' else y0 - side
    return [(x0, y0), (x1, y0), (x1, y1), (x0, y1)]

def make_square(init_pos, side, direction, color=None):
    vertices = square_vertices(init_pos, side, direction)
    x, y = zip(*vertices)
    minx = min(x)
    miny = min(y)
    if not color:
        color = '#'+random.randbytes(3).hex()
    return Rectangle((minx, miny), side, side, color=color, alpha=0.5), vertices

@lru_cache(maxsize=None)
def first_level(number):
    assert number > 2
    half = 180 / number
    x1 = (16 / (1 + tan(half) ** 2)) ** .5
    y1 = (16 - x1 ** 2) ** .5
    k = y1 - tan(half - 90) * x1
    x2 = -k / tan(half - 90)
    r = ((x1 - x2) ** 2 + y1 ** 2) ** .5
    return (x2, r)

@lru_cache(maxsize=None)
def next_level(distance, radius):
    base = (distance - radius) / (distance + radius)
    new_radius = radius * base
    new_distance = distance - radius - new_radius
    return (new_distance, new_radius)

@lru_cache(maxsize=None)
def fade_color(color, strength):
    assert re.match(r'^#[0-9a-f]{6}$', color)
    assert 0 < strength <= 1
    color = int(color[1:], 16)
    r, g, b = [round(((color>>i)%256)*strength) for i in (16, 8, 0)]
    return f'#{r:02x}{g:02x}{b:02x}'

def fibonacci_spiral(n, c=12, bsquare=True, barc=True, reverse=False):
    assert barc or bsquare
    colors = [spectrum_position(int(1530*i/c)) for i in range(c)]
    sequence = fibonacci(n+1)
    init_pos = (0, 0)
    vertexes = []
    shapes = []
    for i, n in enumerate(sequence[1:]):
        direction = DIRECTIONS[i % 4] if not reverse else DIRECTIONS[::-1][i % 4]
        angle = ANGLES[i % 4] if not reverse else CCW_ANGLES[i % 4]
        square, vertices = make_square(init_pos, n, direction, colors[i%c])
        center = vertices[[1, 3][i % 2]]
        a, b = angle
        arc = Arc(center, 2*n, 2*n, theta1=a, theta2=b, edgecolor='k', lw=3.0)
        if bsquare: shapes.append(square)
        if barc: shapes.append(arc)
        init_pos = vertices[-2]
        vertexes.extend(vertices)
    
    return shapes, vertexes

def fibonacci_arcs(n, degree, color='k' ,reverse=False, alpha=0.75, lw=2):
    sequence = fibonacci(n)
    init_pos = (0, 0)
    vertexes = []
    shapes = []
    for i, n in enumerate(sequence[1:]):
        direction = DIRECTIONS[i % 4] if not reverse else DIRECTIONS[::-1][i % 4]
        angle = ANGLES[i % 4] if not reverse else CCW_ANGLES[i % 4]
        vertices = square_vertices(init_pos, n, direction)
        rotated_vertices = [rotate(pos, degree) for pos in vertices]
        index = [1, 3][i % 2]
        center = rotated_vertices[index]
        a, b = [(i + degree) % 360 for i in angle]
        arc = Arc(center, 2*n, 2*n, theta1=a, theta2=b, edgecolor=color, lw=lw, alpha=alpha)
        init_pos = vertices[-2]
        shapes.append(arc)
        vertexes.extend(rotated_vertices)
    
    return shapes, vertexes

################################################################
# plotting functions

def plot_hailstone(n, real=False):
    assert n > 1
    fig = plt.figure(frameon=False)
    ax = fig.add_subplot(111)
    ax.set_axis_off()
    numbers = list(range(-n, 0)) + list(range(1, n))
    for i in numbers:
        ax.plot(hailstone(i))
    
    fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
    if real:
        plt.axis('scaled')
    plt.box(False)
    plt.set_cmap('rainbow')
    plt.show()

def plot_spiral(n, c, barc=True, bsquare=True, reverse=False):
    shapes, vertexes = fibonacci_spiral(n, c, barc=barc, bsquare=bsquare, reverse=reverse)
    x, y = zip(*vertexes)
    fig = plt.figure(frameon=False)
    ax = fig.add_subplot(111)
    ax.set_axis_off()
    for shape in shapes:
        ax.add_patch(shape)
    
    fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
    plt.axis('scaled')
    plt.box(False)
    ax = plt.gca()
    ax.set_xlim([min(x), max(x)])
    ax.set_ylim([min(y), max(y)])
    plt.show()

def plot_lolipop(n, arms, lw=1):
    arcs = []
    points = []
    u = 360/arms
    for i in range(arms):
        shapes, vertices = fibonacci_arcs(n, i*u, color='k', alpha=1, lw=lw)
        arcs.extend(shapes)
        points.extend(vertices)
    
    x, y = zip(*points)
    a, b = points[-2]
    radius = (a*a+b*b)**.5
    fig = plt.figure(frameon=False)
    ax = fig.add_subplot(111)
    ax.set_axis_off()
    for arc in arcs:
        ax.add_patch(arc)
    
    ax.add_patch(Arc((0, 0), 2*radius, 2*radius, theta1=0, theta2=360, edgecolor='k', lw=lw))
    fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
    plt.axis('scaled')
    plt.box(False)
    ax = plt.gca()
    ax.set_xlim([min(x), max(x)])
    ax.set_ylim([min(y), max(y)])
    plt.show()

def plot_sunflower(n, cw, color=True, alpha=0.75):
    assert cw in FIBONACCI and cw <= 144
    i = FIBONACCI.index(cw)
    ccw = FIBONACCI[i+1]
    arcs = []
    points = []
    cwu = 360/cw
    ccwu = 360/ccw
    color_cwu = 1530/cw
    color_ccwu = 1530/ccw
    for i in range(cw):
        color1 = spectrum_position(round(i*color_cwu)) if color else 'k'
        shapes, vertices = fibonacci_arcs(n, i*cwu, color=color1, alpha=alpha)
        arcs.extend(shapes)
        points.extend(vertices)
    
    for i in range(ccw):
        color2 = spectrum_position(round(1530-i*color_ccwu)) if color else 'k'
        shapes, vertices = fibonacci_arcs(n, (90+i*ccwu)%360, color=color2, reverse=True, alpha=alpha)
        arcs.extend(shapes)
        points.extend(vertices)
    
    x, y = zip(*points)
    fig = plt.figure(frameon=False)
    ax = fig.add_subplot(111)
    ax.set_axis_off()
    for arc in arcs:
        ax.add_patch(arc)
    
    fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
    plt.axis('scaled')
    plt.box(False)
    ax = plt.gca()
    ax.set_xlim([min(x), max(x)])
    ax.set_ylim([min(y), max(y)])
    plt.show()

def recur_plot_circle(division, levels, fade=0.9):
    assert division > 2
    fig = plt.figure(frameon=False)
    ax = fig.add_subplot(111)
    ax.set_axis_off()
    rotation_diff = 90 / division
    unit_rotation = 360 / division
    color_increment = 1530 / division
    distance, radius = first_level(division)
    max_distance = distance + radius
    for level in range(levels):
        for index in range(division):
            center = rotate((0, distance), level * rotation_diff + index * unit_rotation)
            color = spectrum_position(round(index * color_increment))
            color = fade_color(color, fade**level)
            ax.add_patch(Circle(center, radius, color=color))
        distance, radius = next_level(distance, radius)
    
    fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
    plt.axis('scaled')
    plt.box(False)
    ax = plt.gca()
    ax.set_xlim([-max_distance, max_distance])
    ax.set_ylim([-max_distance, max_distance])
    plt.show()

def plot_star(number):
    assert number > 2
    fig = plt.figure(frameon=False)
    ax = fig.add_subplot(111)
    ax.set_axis_off()
    unit_rotation = 360 / number
    points = []
    for i in range(number):
        points.append(rotate((0, 4), unit_rotation*i))
    
    for a, b in combinations(points, r=2):
        x, y = zip(*[a, b])
        ax.plot(x, y)
    
    fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
    plt.axis('scaled')
    plt.box(False)
    ax = plt.gca()
    ax.set_xlim([-4, 4])
    ax.set_ylim([-4, 4])
    plt.set_cmap('rainbow')
    plt.show()

def plot_tetragram(n):
    assert n > 2 and not (n - 2) % 3
    averages = [dice_average(i) for i in range(2, n)]
    fig = plt.figure(frameon=False)
    ax = fig.add_subplot(111)
    ax.set_axis_off()
    for a, b, c in zip(averages[::3], averages[1::3], averages[2::3]):
        x, y = zip(*[a, b, c])
        ax.triplot(x, y)
        ax.triplot(negate(x), y)
        ax.triplot(negate(x), negate(y))
        ax.triplot(x, negate(y))
        ax.triplot(y, x)
        ax.triplot(negate(y), x)
        ax.triplot(negate(y), negate(x))
        ax.triplot(y, negate(x))
    
    fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
    plt.axis('scaled')
    plt.box(False)
    plt.set_cmap('rainbow')
    plt.show()

if __name__ == '__main__':
    plot_hailstone(128)
    plot_hailstone(12, 1)
    plot_spiral(24, 6)
    plot_spiral(24, 6, reverse=True)
    plot_lolipop(25, 12)
    for i in (2, 3, 5, 8, 13):
        plot_sunflower(24, i)
    
    for i in (3, 6, 12, 24):
        recur_plot_circle(i, round(i/2))
    
    for i in range(4, 13):
        plot_star(i)
    
    for i in range(14, 32, 3):
        plot_tetragram(i)

All examples:

---

How can the script be improved? I am really new to all this stuff.

Answer 1

Thanks for posting this code; all the cool images really caught my eye.

This is just a mini review since I don't have much to offer in the math arena.

DRY

In the spectrum_position function, some of the same numbers are repeated in the elif lines.

if 0 <= n < 255:
    ...
elif 255 <= n < 510:
    ...
elif 510 <= n < 765:

When the if condition is false, we know that n must be greater than or equal to 255. Therefore, the lower bound check is not needed. We can simplify:

elif 255 <= n < 510:

as:

elif n < 510:

The same is true for the remaining lower bound checks. This makes the code a little easier to read, and it may be slightly more efficient performing fewer checks.

In the plot_tetragram function, these two expressions are repeated several times:

negate(x)
negate(y)

Consider assigning each to a variable:

x, y = zip(*[a, b, c])
neg_x = negate(x)
neg_y = negate(y)
ax.triplot(x, y)
ax.triplot(neg_x, y)

Readability

In lines like this:

radius = (a*a+b*b)**.5

it might be easier to understand the code using the sqrt function:

radius = math.sqrt(a*a+b*b)

This is part of the math library, which you are already importing.

Documentation

It would be nice to have a docstring at the top of the code to summarize its purpose, as well as some sprinkled throughout the code for some of the more complex functions.

Magic number

The number 1530 is used in a few places throughout the code. Consider making it a named constant, for example:

NUM_COLORS = 6*255