7-spheres-3d

Python scripts to generate 3D models of some complex curves

Diffstats

4 files changed, 192 insertions, 0 deletions

diff --git /dev/null b/.gitignore
@@ -0,0 +1,3 @@
+__pycache__
+site/
+models.blend

diff --git /dev/null b/main.py
@@ -0,0 +1,89 @@
+import numpy as np
+from model import Model, normalized_p, riemman_sphere_p
+
+def s(k: int, n: int) -> np.complex128:
+    return np.exp(2 * np.pi * 1j * k / n)
+
+EPS = 0.05
+XI_MAX = 15
+
+def brieskorn_surface(n: int,
+                      eps: float = EPS, xi_max: float = XI_MAX) -> Model:
+    """
+    Returns a model of the Brieskorn surface z_1^3 + z_2^(6 n - 1) = 1
+    """
+
+    if not (1 <= n <= 11):
+        raise ValueError(f"Parameter \"n\" should be in between 1 and 11: n = {n}")
+
+    m = Model(normalized_p)
+    n_ = 6*n - 1
+
+    thetas = np.arange(0, np.pi / 2 + eps, eps)
+    xis = np.arange(-xi_max, xi_max + eps, eps)
+
+    for k1 in range(0, 3):
+        for k2 in range(0, n_):
+            patch = []
+
+            for theta in thetas:
+                row = []
+                for xi in xis:
+                    u1 = (np.exp(xi + 1j * theta)+np.exp(-xi - 1j * theta))/2
+                    z1 = s(k1, 3) * (u1**(2/3))
+                    u2 = (np.exp(xi + 1j * theta)-np.exp(-xi - 1j * theta))/2j
+                    z2 = s(k2, n_) * (u2**(2/n_))
+
+                    row.append(np.array([z1, z2]))
+
+                patch.append(row)
+
+            m.append_patch(patch)
+
+    return m
+
+def fermat_surface(n: int, eps: float = EPS, xi_max: float = XI_MAX) -> Model:
+    """
+    Returns a model of the Fermat surface z_1^n + z_2^n = 1
+    """
+
+    if n <= 0:
+        raise ValueError(f"Parameter \"n\" should be positive: n = {n}")
+
+    m = Model(normalized_p)
+
+    thetas = np.arange(0, np.pi / 2 + eps, eps)
+    xis = np.arange(-xi_max, xi_max + eps, eps)
+
+    for k1 in range(0, n):
+        for k2 in range(0, n):
+            patch = []
+
+            for theta in thetas:
+                row = []
+                for xi in xis:
+                    u1 = (np.exp(xi + 1j * theta)+np.exp(-xi - 1j * theta))/2
+                    z1 = s(k1, n) * (u1**(2/n))
+                    u2 = (np.exp(xi + 1j * theta)-np.exp(-xi - 1j * theta))/2j
+                    z2 = s(k2, n) * (u2**(2/n))
+
+                    row.append(np.array([z1, z2]))
+
+                patch.append(row)
+
+            m.append_patch(patch)
+
+    return m
+
+def main():
+    # print("Creating Brieskorn model for n = 1...")
+    # m = brieskorn_surface(1, eps=0.01)
+    # m.write_obj("models/brieskorn1.obj")
+    # print("Done!")
+
+    print("Creating Brieskorn model for n = 6...")
+    m = brieskorn_surface(6, eps=0.01)
+    m.write_obj("models/brieskorn6.obj")
+    print("Done!")
+
+if __name__ == "__main__": main()

diff --git /dev/null b/model.py
@@ -0,0 +1,100 @@
+import numpy.linalg as la
+import numpy as np
+from typing import List, Tuple, Iterator, Callable
+
+class Model:
+    def __init__(self, p: Callable[[np.complex128, np.complex128], np.array]):
+        self.patches = []
+        self.bounding_x = (0, 0)
+        self.bounding_y = (0, 0)
+        self.bounding_z = (0, 0)
+        self.projection = p
+
+    def append_patch(self, patch: List[List[np.array]]):
+        """
+        Appends a 2D array of points in ℂ² to the model
+        and updates the bounding box
+        """
+        x_min, x_max = self.bounding_x
+        y_min, y_max = self.bounding_y
+        z_min, z_max = self.bounding_z
+
+        patch_real = [[self.projection(z1, z2) for z1, z2 in row]
+                      for row in patch]
+
+        for row in patch_real:
+            for x, y, z in row:
+                x_min = min(x, x_min)
+                x_max = max(x, x_max)
+                y_min = min(y, y_min)
+                y_max = max(y, y_max)
+                z_min = min(z, z_min)
+                z_max = max(z, z_max)
+
+        self.patches.append(patch_real)
+        self.bounding_x = x_min, x_max 
+        self.bounding_y = y_min, y_max 
+        self.bounding_z = z_min, z_max 
+
+    def vertices(self) -> Iterator[Tuple[np.float64]]:
+        """
+        Iterate through the vertices of the mesh in normalized position:
+        vertices are aligned to the center using the model's bounding box
+        """
+        x_min, x_max = self.bounding_x
+        y_min, y_max = self.bounding_y
+        z_min, z_max = self.bounding_z
+
+        x_off = -(x_min + x_max)/2
+        y_off = -(y_min + y_max)/2
+        z_off = -(z_min + z_max)/2
+
+        return ((x + x_off, y + y_off, z + z_off)
+                for patch in self.patches for row in patch for x, y, z in row)
+
+    def write_obj(self, path: str):
+        with open(path, "w") as f:
+            f.write("# vertices\n")
+
+            for x, y, z in self.vertices():
+                f.write(f"v {x:.3f} {y:.3f} {z:.3f}\n")
+
+            f.write("\n# faces\n")
+
+            v_i = 1 # OBJ indices start at 1
+            for patch in self.patches:
+                rows = len(patch)
+                cols = len(patch[0])
+
+                for theta_i in range(rows - 1):
+                    for xi_i in range(cols - 1):
+                        v1 = v_i + xi_i + theta_i * cols
+                        v2 = v_i + xi_i + (theta_i + 1) * cols
+                        v3 = v_i + (xi_i + 1) + theta_i * cols
+                        v4 = v_i + (xi_i + 1) + (theta_i + 1) * cols
+
+                        f.write(f"f {v1} {v2} {v3}\n")
+                        f.write(f"f {v2} {v4} {v3}\n")
+
+                v_i += rows * cols
+
+def normalized_p(z1: np.complex128, z2: np.complex128) -> np.array:
+    "Projection ℝ⁴ → 𝕊³ + inclusion ℝ³ → 𝕊³"
+    v = np.array([z1, z2])
+    v = v/la.norm(v)
+
+    v_1, v_2 = v
+    x, y = v_1.real, v_1.imag
+    z, w = v_2.real, v_2.imag
+
+    return np.array([x, y, z]) / (1 - w)
+
+def riemman_sphere_p(z1: np.complex128, z2: np.complex128) -> np.array:
+    "Riemman-sphere transform + projection on the first 3 coordinates"
+    v = np.array([z1.real, z1.imag, z2.real, z2.imag])
+    n_sqrd = la.norm(v)**2
+
+    u1, u2, _, _ = v / (1 + n_sqrd/4)
+    u0 = 2 * (n_sqrd/4)/(1 + n_sqrd/4)
+
+    return np.array([u0, u1, u2])

diff --git /dev/null b/references/Hanson - A Construction for Computer Visualization of Certain Complex Curves.pdf
Binary files differ