diff --git a/main.py b/main.py
index 42fecfd..b76615c 100644
--- a/main.py
+++ b/main.py
@@ -1,13 +1,11 @@
 import io
-import os
-import matplotlib
+from perlin import CoordsGenerator
 from typing import Literal
 import itertools
 import logging
 import random
 
 import matplotlib.pyplot as plt
-import matplotlib.font_manager as font_manager
 import numpy as np
 import uvicorn
 from fastapi import FastAPI
@@ -90,12 +88,13 @@ marker_subs = {
 @app.get("/speckles/")
 def make_wallpaper(
     speckle_colours: str,
-    density: float | None = 0.12,
-    size: float | None = 3,
+    density: float = 0.12,
+    size: float = 3,
     fileformat: str = "svg",
-    orientation: str | None = "landscape",
+    orientation: str = "landscape",
     local: bool = False,
-    markers: str | None = ".",
+    markers: str = ".",
+    perlin: bool = True,
 ):
     if fileformat not in MEDIA_TYPES:
         return
@@ -126,15 +125,25 @@ def make_wallpaper(
     ax.set_yticks([])
     ax.margins(0, 0)
 
+    if perlin:
+        gen = CoordsGenerator(y, x)
+
     for color, marker, size in itertools.product(
         speckle_colours,
         markers,
         np.logspace(0, size, 10, base=np.exp(2)),
     ):
         marker = marker_subs.get(marker, marker)
+        if perlin:
+            x_coords, y_coords = gen.pick(speckles_per_colour)
+        else:
+            x_coords, y_coords = (
+                [random.random() * x / 8 for _ in range(speckles_per_colour)],
+                [random.random() * y / 8 for _ in range(speckles_per_colour)],
+            )
         ax.scatter(
-            [random.random() * x / 8 for _ in range(speckles_per_colour)],
-            [random.random() * y / 8 for _ in range(speckles_per_colour)],
+            x_coords,
+            y_coords,
             c=color,
             s=size,
             marker=marker,
diff --git a/perlin.py b/perlin.py
new file mode 100644
index 0000000..4d47bd9
--- /dev/null
+++ b/perlin.py
@@ -0,0 +1,84 @@
+import numpy as np
+from scipy.special import softmax
+
+
+def interpolant(t):
+    return t * t * t * (t * (t * 6 - 15) + 10)
+
+
+def generate_perlin_noise_2d(
+    shape, res, tileable=(False, False), interpolant=interpolant
+):
+    """Generate a 2D numpy array of perlin noise.
+
+    Args:
+        shape: The shape of the generated array (tuple of two ints).
+            This must be a multple of res.
+        res: The number of periods of noise to generate along each
+            axis (tuple of two ints). Note shape must be a multiple of
+            res.
+        tileable: If the noise should be tileable along each axis
+            (tuple of two bools). Defaults to (False, False).
+        interpolant: The interpolation function, defaults to
+            t*t*t*(t*(t*6 - 15) + 10).
+
+    Returns:
+        A numpy array of shape shape with the generated noise.
+
+    Raises:
+        ValueError: If shape is not a multiple of res.
+    """
+    delta = (res[0] / shape[0], res[1] / shape[1])
+    d = (shape[0] // res[0], shape[1] // res[1])
+    grid = np.mgrid[0 : res[0] : delta[0], 0 : res[1] : delta[1]].transpose(1, 2, 0) % 1
+    # Gradients
+    angles = 2 * np.pi * np.random.rand(res[0] + 1, res[1] + 1)
+    gradients = np.dstack((np.cos(angles), np.sin(angles)))
+    if tileable[0]:
+        gradients[-1, :] = gradients[0, :]
+    if tileable[1]:
+        gradients[:, -1] = gradients[:, 0]
+    gradients = gradients.repeat(d[0], 0).repeat(d[1], 1)
+    g00 = gradients[: -d[0], : -d[1]]
+    g10 = gradients[d[0] :, : -d[1]]
+    g01 = gradients[: -d[0], d[1] :]
+    g11 = gradients[d[0] :, d[1] :]
+    # Ramps
+    n00 = np.sum(np.dstack((grid[:, :, 0], grid[:, :, 1])) * g00, 2)
+    n10 = np.sum(np.dstack((grid[:, :, 0] - 1, grid[:, :, 1])) * g10, 2)
+    n01 = np.sum(np.dstack((grid[:, :, 0], grid[:, :, 1] - 1)) * g01, 2)
+    n11 = np.sum(np.dstack((grid[:, :, 0] - 1, grid[:, :, 1] - 1)) * g11, 2)
+    # Interpolation
+    t = interpolant(grid)
+    n0 = n00 * (1 - t[:, :, 0]) + t[:, :, 0] * n10
+    n1 = n01 * (1 - t[:, :, 0]) + t[:, :, 0] * n11
+    return np.sqrt(2) * ((1 - t[:, :, 1]) * n0 + t[:, :, 1] * n1)
+
+
+class CoordsGenerator:
+    def __init__(self, x: int, y: int, factor: int = 500):
+        print(x, y, x // factor, y // factor)
+        self.noise = generate_perlin_noise_2d((x, y), (x // factor, y // factor))
+        self.noise_distribution = softmax(self.noise, axis=1)
+
+    def pick(self, n):
+        x, y = self.noise.shape
+        x = np.random.choice(x, size=n, replace=False)
+        y = [
+            np.random.choice(y, size=1, p=self.noise_distribution[x_, :], replace=False) for x_ in x
+        ]
+        return x, y
+
+
+if __name__ == "__main__":
+    import matplotlib.pyplot as plt
+
+    # gen = CoordsGenerator(1920, 1080, threshold=0.85)
+    # x, y = gen.pick(1000)
+    # plt.scatter(x, y)
+    factor = 500
+    noise = generate_perlin_noise_2d((1080, 1920), (1080 // factor, 1920 // factor))
+
+    plt.matshow(noise, cmap="bwr")
+    plt.colorbar()
+    plt.show()