Noisy planets

I recently came across this interesting post and I wanted to replicate the whole thing in python within a jupyter notebook.

Let's start by correctly generating random positions within a disk (the original post missed the right way to do it as explained here). We then plot them in polar coordinates.

import numpy as np

npoints = 1000
radius = 10

r = np.sqrt(np.random.rand(npoints)) * radius
theta = np.random.rand(npoints) * 2 * np.pi

import matplotlib.pyplot as plt

fig, ax = plt.subplots(subplot_kw={"projection": "polar"})
ax.plot(theta, r, "o", alpha=0.5);

Given this set of initial point, we will generate random trajectories using Perlin noise. Despite this kind of noise seems not to be coded within numpy or scipy, there are python packages to do it and we will install and use the following one

!pip install perlin-noise

Show code cell output Hide code cell output

Defaulting to user installation because normal site-packages is not writeable
Requirement already satisfied: perlin-noise in /home/garrido/.local/lib/python3.11/site-packages (1.12)

Let's write a function to return a bunch of trajectories cutting each of them when the next point is outside the planet.

from perlin_noise import PerlinNoise


def generate_noisy_planet(npoints=1000, radius=10, step=0.1, seed=None):
    np.random.seed(seed)
    noise = PerlinNoise(octaves=4, seed=seed)

    trajectories = []
    for i in range(npoints):
        r = np.sqrt(np.random.rand()) * radius
        theta = np.random.rand() * 2 * np.pi
        x, y = r * np.cos(theta), r * np.sin(theta)

        trajectory = []
        while np.hypot(x, y) < radius:
            trajectory += [[x, y]]
            n = noise([x / 100, y / 100])
            fx = np.sin if n > 0.5 else np.cos
            fy = np.cos if n > 0.5 else np.sin
            x += fx(n * 2 * np.pi) * step
            y += fy(n * 2 * np.pi) * step

        trajectories += [np.array(trajectory)]
    return trajectories

Finally let's plot them in cartesian coordinates

trajectories = generate_noisy_planet()

fig, ax = plt.subplots(figsize=(8, 8))
ax.axis("equal")
ax.axis("off")
cmap = plt.get_cmap("Greys", len(trajectories))
for i, trajectory in enumerate(trajectories):
    plt.plot(trajectory[:, 0], trajectory[:, 1], color=cmap(i))

And here are some nice examples of noisy planets

seeds = [20, 52, 666]

cmaps = ["Reds", "Greens", "Blues"]
fig, axes = plt.subplots(ncols=3, figsize=(24, 8))

for i, seed in enumerate(seeds):
    axes[i].axis("equal")
    axes[i].axis("off")
    trajectories = generate_noisy_planet(seed=seed)
    cmap = plt.get_cmap(cmaps[i], len(trajectories))
    for j, trajectory in enumerate(trajectories):
        axes[i].plot(*trajectory.T, color=cmap(j))

Just like the original post, we finally create a gif file for the latest evil planet (seed=666)

import os

from tqdm.auto import tqdm

tmp_fig = "/tmp/planets"
os.makedirs(tmp_fig, exist_ok=True)

fig, ax = plt.subplots(figsize=(8, 8))
max_time = np.max([len(trajectory) for trajectory in trajectories])
for t in tqdm(range(0, max_time)):
    ax.axis("equal")
    ax.axis("off")
    for i, trajectory in enumerate(trajectories):
        tlim = min(t, len(trajectory))
        ax.plot(trajectory[:tlim, 0], trajectory[:tlim, 1], color=cmap(i))
    fig.savefig(os.path.join(tmp_fig, f"planet_{t:03d}.png"))
    ax.clear()
plt.close()

and the process conversion using ImageMagick

import subprocess

subprocess.run(
    [
        "convert",
        "-background",
        "white",
        "-alpha",
        "remove",
        "-layers",
        "optimize",
        "-delay",
        "10",
        "-loop",
        "0",
        os.path.join(tmp_fig, "*.png"),
        "./images/animation.gif",
    ]
);

version control - commit 38d018f - 2023-11-01