init

Pipfile

@@ -0,0 +1,14 @@
+[[source]]
+url = "https://pypi.org/simple"
+verify_ssl = true
+name = "pypi"
+
+[packages]
+numpy = "*"
+pytest = "*"
+matplotlib = "*"
+
+[dev-packages]
+
+[requires]
+python_version = "3.7"

Pipfile.lock

@@ -0,0 +1,174 @@
+{
+    "_meta": {
+        "hash": {
+            "sha256": "38985b134d3f8d2a44b9c4d413de1240fd302ad89956c4c8265719c3445ee91b"
+        },
+        "pipfile-spec": 6,
+        "requires": {
+            "python_version": "3.7"
+        },
+        "sources": [
+            {
+                "name": "pypi",
+                "url": "https://pypi.org/simple",
+                "verify_ssl": true
+            }
+        ]
+    },
+    "default": {
+        "atomicwrites": {
+            "hashes": [
+                "sha256:03472c30eb2c5d1ba9227e4c2ca66ab8287fbfbbda3888aa93dc2e28fc6811b4",
+                "sha256:75a9445bac02d8d058d5e1fe689654ba5a6556a1dfd8ce6ec55a0ed79866cfa6"
+            ],
+            "version": "==1.3.0"
+        },
+        "attrs": {
+            "hashes": [
+                "sha256:69c0dbf2ed392de1cb5ec704444b08a5ef81680a61cb899dc08127123af36a79",
+                "sha256:f0b870f674851ecbfbbbd364d6b5cbdff9dcedbc7f3f5e18a6891057f21fe399"
+            ],
+            "version": "==19.1.0"
+        },
+        "cycler": {
+            "hashes": [
+                "sha256:1d8a5ae1ff6c5cf9b93e8811e581232ad8920aeec647c37316ceac982b08cb2d",
+                "sha256:cd7b2d1018258d7247a71425e9f26463dfb444d411c39569972f4ce586b0c9d8"
+            ],
+            "version": "==0.10.0"
+        },
+        "kiwisolver": {
+            "hashes": [
+                "sha256:0ee4ed8b3ae8f5f712b0aa9ebd2858b5b232f1b9a96b0943dceb34df2a223bc3",
+                "sha256:0f7f532f3c94e99545a29f4c3f05637f4d2713e7fd91b4dd8abfc18340b86cd5",
+                "sha256:1a078f5dd7e99317098f0e0d490257fd0349d79363e8c923d5bb76428f318421",
+                "sha256:1aa0b55a0eb1bd3fa82e704f44fb8f16e26702af1a073cc5030eea399e617b56",
+                "sha256:2874060b91e131ceeff00574b7c2140749c9355817a4ed498e82a4ffa308ecbc",
+                "sha256:379d97783ba8d2934d52221c833407f20ca287b36d949b4bba6c75274bcf6363",
+                "sha256:3b791ddf2aefc56382aadc26ea5b352e86a2921e4e85c31c1f770f527eb06ce4",
+                "sha256:4329008a167fac233e398e8a600d1b91539dc33c5a3eadee84c0d4b04d4494fa",
+                "sha256:45813e0873bbb679334a161b28cb9606d9665e70561fd6caa8863e279b5e464b",
+                "sha256:53a5b27e6b5717bdc0125338a822605084054c80f382051fb945d2c0e6899a20",
+                "sha256:574f24b9805cb1c72d02b9f7749aa0cc0b81aa82571be5201aa1453190390ae5",
+                "sha256:66f82819ff47fa67a11540da96966fb9245504b7f496034f534b81cacf333861",
+                "sha256:79e5fe3ccd5144ae80777e12973027bd2f4f5e3ae8eb286cabe787bed9780138",
+                "sha256:83410258eb886f3456714eea4d4304db3a1fc8624623fc3f38a487ab36c0f653",
+                "sha256:8b6a7b596ce1d2a6d93c3562f1178ebd3b7bb445b3b0dd33b09f9255e312a965",
+                "sha256:9576cb63897fbfa69df60f994082c3f4b8e6adb49cccb60efb2a80a208e6f996",
+                "sha256:95a25d9f3449046ecbe9065be8f8380c03c56081bc5d41fe0fb964aaa30b2195",
+                "sha256:a424f048bebc4476620e77f3e4d1f282920cef9bc376ba16d0b8fe97eec87cde",
+                "sha256:aaec1cfd94f4f3e9a25e144d5b0ed1eb8a9596ec36d7318a504d813412563a85",
+                "sha256:acb673eecbae089ea3be3dcf75bfe45fc8d4dcdc951e27d8691887963cf421c7",
+                "sha256:b15bc8d2c2848a4a7c04f76c9b3dc3561e95d4dabc6b4f24bfabe5fd81a0b14f",
+                "sha256:b1c240d565e977d80c0083404c01e4d59c5772c977fae2c483f100567f50847b",
+                "sha256:c595693de998461bcd49b8d20568c8870b3209b8ea323b2a7b0ea86d85864694",
+                "sha256:ce3be5d520b4d2c3e5eeb4cd2ef62b9b9ab8ac6b6fedbaa0e39cdb6f50644278",
+                "sha256:e0f910f84b35c36a3513b96d816e6442ae138862257ae18a0019d2fc67b041dc",
+                "sha256:ea36e19ac0a483eea239320aef0bd40702404ff8c7e42179a2d9d36c5afcb55c",
+                "sha256:efabbcd4f406b532206b8801058c8bab9e79645b9880329253ae3322b7b02cd5",
+                "sha256:f923406e6b32c86309261b8195e24e18b6a8801df0cfc7814ac44017bfcb3939"
+            ],
+            "version": "==1.0.1"
+        },
+        "matplotlib": {
+            "hashes": [
+                "sha256:1ae6549976b6ceb6ee426272a28c0fc9715b3e3669694d560c8f661c5b39e2c5",
+                "sha256:4d4250bf508dd07cca3b43888097f873cadb66eec6ac63dbbfb798798ec07af2",
+                "sha256:53af2e01d7f1700ed2b64a9091bc865360c9c4032f625451c4589a826854c787",
+                "sha256:63e498067d32d627111cd1162cae1621f1221f9d4c6a9745dd7233f29de581b6",
+                "sha256:7169a34971e398dd58e87e173f97366fd88a3fa80852704530433eb224a8ca57",
+                "sha256:91c54d6bb9eeaaff965656c5ea6cbdcbf780bad8462ac99b30b451548194746f",
+                "sha256:aeef177647bb3fccfe09065481989d7dfc5ac59e9367d6a00a3481062cf651e4",
+                "sha256:cf8ae10559a78aee0409ede1e9d4fda03895433eeafe609dd9ed67e45f552db0",
+                "sha256:d51d0889d1c4d51c51a9822265c0494ea3e70a52bdd88358e0863daca46fa23a",
+                "sha256:de5ccd3500247f85fe4f9fad90f80a8bd397e4f110a4c33fabf95f07403e8372",
+                "sha256:e1d33589e32f482d0a7d1957bf473d43341115d40d33f578dad44432e47df7b7",
+                "sha256:e8d1939262aa6b36d0c51f50a50a43a04b9618d20db31e6c0192b1463067aeef",
+                "sha256:e918d51b1fda82a65fdf52d2f3914b2246481cc2a9cd10e223e6be6078916ff3"
+            ],
+            "index": "pypi",
+            "version": "==3.0.3"
+        },
+        "more-itertools": {
+            "hashes": [
+                "sha256:0125e8f60e9e031347105eb1682cef932f5e97d7b9a1a28d9bf00c22a5daef40",
+                "sha256:590044e3942351a1bdb1de960b739ff4ce277960f2425ad4509446dbace8d9d1"
+            ],
+            "markers": "python_version > '2.7'",
+            "version": "==6.0.0"
+        },
+        "numpy": {
+            "hashes": [
+                "sha256:1980f8d84548d74921685f68096911585fee393975f53797614b34d4f409b6da",
+                "sha256:22752cd809272671b273bb86df0f505f505a12368a3a5fc0aa811c7ece4dfd5c",
+                "sha256:23cc40313036cffd5d1873ef3ce2e949bdee0646c5d6f375bf7ee4f368db2511",
+                "sha256:2b0b118ff547fecabc247a2668f48f48b3b1f7d63676ebc5be7352a5fd9e85a5",
+                "sha256:3a0bd1edf64f6a911427b608a894111f9fcdb25284f724016f34a84c9a3a6ea9",
+                "sha256:3f25f6c7b0d000017e5ac55977a3999b0b1a74491eacb3c1aa716f0e01f6dcd1",
+                "sha256:4061c79ac2230594a7419151028e808239450e676c39e58302ad296232e3c2e8",
+                "sha256:560ceaa24f971ab37dede7ba030fc5d8fa173305d94365f814d9523ffd5d5916",
+                "sha256:62be044cd58da2a947b7e7b2252a10b42920df9520fc3d39f5c4c70d5460b8ba",
+                "sha256:6c692e3879dde0b67a9dc78f9bfb6f61c666b4562fd8619632d7043fb5b691b0",
+                "sha256:6f65e37b5a331df950ef6ff03bd4136b3c0bbcf44d4b8e99135d68a537711b5a",
+                "sha256:7a78cc4ddb253a55971115f8320a7ce28fd23a065fc33166d601f51760eecfa9",
+                "sha256:80a41edf64a3626e729a62df7dd278474fc1726836552b67a8c6396fd7e86760",
+                "sha256:893f4d75255f25a7b8516feb5766c6b63c54780323b9bd4bc51cdd7efc943c73",
+                "sha256:972ea92f9c1b54cc1c1a3d8508e326c0114aaf0f34996772a30f3f52b73b942f",
+                "sha256:9f1d4865436f794accdabadc57a8395bd3faa755449b4f65b88b7df65ae05f89",
+                "sha256:9f4cd7832b35e736b739be03b55875706c8c3e5fe334a06210f1a61e5c2c8ca5",
+                "sha256:adab43bf657488300d3aeeb8030d7f024fcc86e3a9b8848741ea2ea903e56610",
+                "sha256:bd2834d496ba9b1bdda3a6cf3de4dc0d4a0e7be306335940402ec95132ad063d",
+                "sha256:d20c0360940f30003a23c0adae2fe50a0a04f3e48dc05c298493b51fd6280197",
+                "sha256:d3b3ed87061d2314ff3659bb73896e622252da52558f2380f12c421fbdee3d89",
+                "sha256:dc235bf29a406dfda5790d01b998a1c01d7d37f449128c0b1b7d1c89a84fae8b",
+                "sha256:fb3c83554f39f48f3fa3123b9c24aecf681b1c289f9334f8215c1d3c8e2f6e5b"
+            ],
+            "index": "pypi",
+            "version": "==1.16.2"
+        },
+        "pluggy": {
+            "hashes": [
+                "sha256:19ecf9ce9db2fce065a7a0586e07cfb4ac8614fe96edf628a264b1c70116cf8f",
+                "sha256:84d306a647cc805219916e62aab89caa97a33a1dd8c342e87a37f91073cd4746"
+            ],
+            "version": "==0.9.0"
+        },
+        "py": {
+            "hashes": [
+                "sha256:64f65755aee5b381cea27766a3a147c3f15b9b6b9ac88676de66ba2ae36793fa",
+                "sha256:dc639b046a6e2cff5bbe40194ad65936d6ba360b52b3c3fe1d08a82dd50b5e53"
+            ],
+            "version": "==1.8.0"
+        },
+        "pyparsing": {
+            "hashes": [
+                "sha256:66c9268862641abcac4a96ba74506e594c884e3f57690a696d21ad8210ed667a",
+                "sha256:f6c5ef0d7480ad048c054c37632c67fca55299990fff127850181659eea33fc3"
+            ],
+            "version": "==2.3.1"
+        },
+        "pytest": {
+            "hashes": [
+                "sha256:592eaa2c33fae68c7d75aacf042efc9f77b27c08a6224a4f59beab8d9a420523",
+                "sha256:ad3ad5c450284819ecde191a654c09b0ec72257a2c711b9633d677c71c9850c4"
+            ],
+            "index": "pypi",
+            "version": "==4.3.1"
+        },
+        "python-dateutil": {
+            "hashes": [
+                "sha256:7e6584c74aeed623791615e26efd690f29817a27c73085b78e4bad02493df2fb",
+                "sha256:c89805f6f4d64db21ed966fda138f8a5ed7a4fdbc1a8ee329ce1b74e3c74da9e"
+            ],
+            "version": "==2.8.0"
+        },
+        "six": {
+            "hashes": [
+                "sha256:3350809f0555b11f552448330d0b52d5f24c91a322ea4a15ef22629740f3761c",
+                "sha256:d16a0141ec1a18405cd4ce8b4613101da75da0e9a7aec5bdd4fa804d0e0eba73"
+            ],
+            "version": "==1.12.0"
+        }
+    },
+    "develop": {}
+}

src/geometry/affine.py

@@ -0,0 +1,48 @@
+from dataclasses import dataclass
+
+import numpy as np
+
+from geometry.vector import Vector
+
+
+@dataclass
+class AffineTransform:
+  A: np.ndarray
+  A_extended: np.ndarray
+
+  inv_A: np.ndarray
+  inv_A_extended: np.ndarray
+
+  def __init__(self, A: np.ndarray, b: Vector):
+    n = len(A)
+
+    assert A.shape == (n, n)
+    assert b.a.shape == (n, 1)
+
+    self.A = np.array(A, dtype=float)
+    self.A.setflags(write=False)
+
+    self.A_extended = np.zeros((n + 1, n + 1), dtype=float)
+    self.A_extended[:-1, :-1] = self.A
+    self.A_extended[:-1, -1:] = b.a
+    self.A_extended[-1, -1] = 1
+    self.A_extended.setflags(write=False)
+
+    self.inv_A = np.linalg.inv(self.A)
+    self.inv_A_extended = np.linalg.inv(self.A_extended)
+
+  def apply_shift(self, v: Vector) -> Vector:
+    v_extended = np.vstack((v.a, np.eye(1)))
+    result_extended = self.A_extended @ v_extended
+    return Vector(result_extended[:-1, :])
+
+  def apply_no_shift(self, v: Vector) -> Vector:
+    return Vector(self.A @ v.a)
+
+  def invert_shift(self, v: Vector) -> Vector:
+    v_extended = np.vstack((v.a, np.eye(1)))
+    result_extended = self.inv_A_extended @ v_extended
+    return Vector(result_extended[:-1, :])
+
+  def invert_no_shift(self, v: Vector) -> Vector:
+    return Vector(self.inv_A @ v.a)

src/geometry/ellipse.py

@@ -0,0 +1,167 @@
+from dataclasses import dataclass
+from typing import Optional
+
+import numpy as np
+
+from geometry.ray import Ray
+from geometry.surface import Surface
+from geometry.line_drawers import draw_line_3d
+from geometry.vector import Vector
+from geometry.rotations import rotate2d, rotate3d, rotate_rodrigues
+from util import double_equality, minpos, solve_quadratic
+
+
+@dataclass
+class Ellipse(Surface):
+  refraction_index: float
+  f1: Vector
+  f2: Vector
+  distance: float
+
+  def __post_init__(self):
+    if abs(self.f1 - self.f2) > self.distance:
+      raise ValueError()
+
+  def distance_to(self, v):
+    return abs(self.f1 - v) + abs(self.f2 - v)
+
+  def __and__(self, v: Vector) -> bool:
+    return double_equality(self.distance_to(v), self.distance)
+
+  def __contains__(self, v: Vector) -> bool:
+    return self & v or self.distance_to(v) < self.distance
+
+  def normal_at(self, v: Vector) -> Vector:
+    d1 = (v - self.f1).normalize()
+    d2 = (v - self.f2).normalize()
+    return (d1 + d2).normalize()
+
+  def intersect(self, ray: Ray) -> Optional[float]:
+    D = self.distance
+    b1 = 2 * (ray.e @ (ray.r0 - self.f1))
+    b2 = 2 * (ray.e @ (ray.r0 - self.f2))
+    c1 = (ray.r0 - self.f1) @ (ray.r0 - self.f1)
+    c2 = (ray.r0 - self.f2) @ (ray.r0 - self.f2)
+
+    a = 4 * D ** 2 - (b1 - b2) ** 2
+    b = 2 * D ** 2 * (b1 + b2) - 2 * b1 * c1 + 2 * b1 * c2 + 2 * b2 * c1 - 2 * b2 * c2
+    c = 2 * D ** 2 * (c1 + c2) - (c1 - c2) ** 2 - D ** 4
+
+    return minpos(solve_quadratic(a, b, c))
+
+  def draw2d(self, plt, **kwargs):
+    r0 = (self.f1 + self.f2) / 2
+    e = Vector.fromargs(1, 0)
+
+    xs = []
+    ys = []
+
+    for alpha in np.linspace(0, 2 * np.pi, 200):
+      r = Ray(r0=r0, e=rotate2d(e, alpha))
+      p = r.prolong(self.intersect(r))
+      xs.append(p.a[0, 0])
+      ys.append(p.a[1, 0])
+
+      # n = self.normal_at(p)
+      # p2 = p + n
+      # plt.plot(
+      #   [p.a[0, 0], p2.a[0, 0]],
+      #   [p.a[1, 0], p2.a[1, 0]],
+      # )
+
+    plt.plot(xs, ys, **kwargs)
+
+  def draw3d_old(self, axis, **kwargs):
+    r0 = (self.f1 + self.f2) / 2
+    e = Vector.fromargs(0, 0, 1)
+
+    N1 = 20
+    N2 = 100
+
+    for alpha in np.linspace(0, 2 * np.pi, N1):
+      xs = []
+      ys = []
+      zs = []
+
+      for beta in np.linspace(0, 2 * np.pi, N2):
+        d = Vector.fromargs(
+          np.sin(alpha) * np.cos(beta),
+          np.sin(alpha) * np.sin(beta),
+          np.cos(alpha),
+        )
+
+        r = Ray(r0=r0, e=d)
+        p = r.prolong(self.intersect(r))
+
+        xs.append(p.a[0, 0])
+        ys.append(p.a[1, 0])
+        zs.append(p.a[2, 0])
+
+      axis.plot(xs, ys, zs, color='blue', **kwargs)
+
+      xs = []
+      ys = []
+      zs = []
+
+      for beta in np.linspace(0, 2 * np.pi, N2):
+        d = Vector.fromargs(
+          np.sin(beta) * np.cos(alpha),
+          np.sin(beta) * np.sin(alpha),
+          np.cos(beta),
+        )
+
+        r = Ray(r0=r0, e=d)
+        p = r.prolong(self.intersect(r))
+
+        xs.append(p.a[0, 0])
+        ys.append(p.a[1, 0])
+        zs.append(p.a[2, 0])
+
+      axis.plot(xs, ys, zs, color='blue', **kwargs)
+
+  def draw3d(self, axis, **kwargs):
+    r0 = (self.f1 + self.f2) / 2
+    a1 = (self.f2 - r0).normalize()
+
+    a2_candidates = [
+      a1 | Vector.fromargs(1, 0, 0),
+      a1 | Vector.fromargs(0, 1, 0),
+      a1 | Vector.fromargs(0, 0, 1),
+    ]
+    a2 = a2_candidates[np.argmax(np.abs(a2_candidates))]
+    a2 = a2.normalize()
+    a3 = (a1 | a2).normalize()
+
+    draw_line_3d(axis, r0, r0 + a1, lw=5, zorder=0, color='red')
+    draw_line_3d(axis, r0, r0 + a2, lw=5, zorder=0, color='green')
+    draw_line_3d(axis, r0, r0 + a3, lw=5, zorder=0, color='blue')
+
+    N1 = 8
+    N2 = 100
+
+    for lat in np.linspace(0, np.pi, N1):
+      vs = []
+      e = rotate_rodrigues(a1, a2, lat)
+
+      for long in np.linspace(0, 2 * np.pi, N2):
+        r = Ray(r0=r0, e=rotate_rodrigues(e, a1, long))
+        p = r.prolong(self.intersect(r))
+
+        vs.append(p.a)
+
+      V = np.hstack(vs)
+      axis.plot(V[0, :], V[1, :], V[2, :], color='navy', **kwargs)
+
+    for long in np.linspace(0, np.pi, N1):
+      vs = []
+      e = rotate_rodrigues(a2, a1, long)
+      a3_e = rotate_rodrigues(a3, a1, long)
+
+      for lat in np.linspace(0, 2 * np.pi, N2):
+        r = Ray(r0=r0, e=rotate_rodrigues(e, a3_e, lat))
+        p = r.prolong(self.intersect(r))
+
+        vs.append(p.a)
+
+      V = np.hstack(vs)
+      axis.plot(V[0, :], V[1, :], V[2, :], color='indigo', **kwargs)

src/geometry/line_drawers.py

@@ -0,0 +1,55 @@
+from dataclasses import dataclass
+from typing import List
+
+import numpy as np
+from matplotlib import pyplot as plt
+
+from geometry.affine import AffineTransform
+from geometry.vector import Vector
+
+
+def draw_line_2d(axis, v1: Vector, v2: Vector, **kwargs):
+  axis.plot(
+    [v1.a[0, 0], v2.a[0, 0]],
+    [v1.a[1, 0], v2.a[1, 0]],
+    **kwargs)
+
+
+def draw_line_3d(axis, v1: Vector, v2: Vector, **kwargs):
+  axis.plot(
+    [v1.a[0, 0], v2.a[0, 0]],
+    [v1.a[1, 0], v2.a[1, 0]],
+    [v1.a[2, 0], v2.a[2, 0]],
+    **kwargs)
+
+
+class Painter:
+  axis: plt.Axes
+
+  def __init__(self, axis):
+    self.axis = axis
+
+  def plot3d(self, vs: List[Vector], **kwargs):
+    V = np.hstack([v.a for v in vs])
+
+    self.axis.plot(V[0, :], V[1, :], V[2, :], **kwargs)
+
+  def line3d(self, a: Vector, b: Vector, **kwargs):
+    self.plot3d([a, b], **kwargs)
+
+
+class TransformedPainter:
+  painter: Painter
+  transform: AffineTransform
+
+  def __init__(self, painter, transform):
+    self.painter = painter
+    self.transform = transform
+
+  def plot3d(self, vs: List[Vector], **kwargs):
+    vs_ = list(map(self.transform.invert_shift, vs))
+
+    self.painter.plot3d(vs_, **kwargs)
+
+  def line3d(self, a: Vector, b: Vector, **kwargs):
+    self.plot3d([a, b], **kwargs)

src/geometry/plane.py

@@ -0,0 +1,35 @@
+from dataclasses import dataclass
+from typing import Optional
+
+from ray import Ray
+from surface import Surface
+from util import minpos
+from vector import Vector
+
+
+@dataclass
+class Plane(Surface):
+  refraction_index: float
+  r0: Vector
+  normal: Vector
+
+  def __post_init__(self):
+    self.normal = self.normal.normalize()
+
+  def __and__(self, v: Vector) -> bool:
+    return self.normal % (self.r0 - v)
+
+  def __contains__(self, v: Vector) -> bool:
+    proj = (self.r0 - v) @ self.normal
+    return self & v or proj > 0
+
+  def normal_at(self, r: Vector) -> Vector:
+    return self.normal
+
+  def intersect(self, ray: Ray) -> Optional[float]:
+    if self.normal % ray.e:
+      return None
+
+    t = (self.normal @ (self.r0 - ray.r0)) / (self.normal @ ray.e)
+
+    return minpos([t])

src/geometry/ray.py

@@ -0,0 +1,15 @@
+from dataclasses import dataclass
+
+from geometry.vector import Vector
+
+
+@dataclass
+class Ray:
+  r0: Vector
+  e: Vector
+
+  def __post_init__(self):
+    self.e = self.e.normalize()
+
+  def prolong(self, t: float) -> Vector:
+    return self.r0 + self.e * t

src/geometry/reflections.py

@@ -0,0 +1,20 @@
+import numpy as np
+
+from geometry.vector import Vector
+
+
+def reflect(direction: Vector, normal: Vector) -> Vector:
+  reflected = direction - 2 * (direction @ normal) * normal
+  return reflected.normalize()
+
+
+def refract(direction: Vector, normal: Vector, n1: float, n2: float) -> Vector:
+  if n1 > n2:
+    return refract(direction, -1 * normal, n2, n1)
+
+  proj = normal @ (n1 * direction)
+  k = (n2 ** 2 - n1 ** 2) / proj ** 2 + 1
+  refracted = (1 / n2) * (n1 * direction + normal * proj
+                          * (np.sqrt(k) - 1))
+
+  return refracted.normalize()

src/geometry/rotations.py

@@ -0,0 +1,45 @@
+import numpy as np
+
+from geometry.vector import Vector
+
+
+def rotate2d(v: Vector, alpha: float) -> Vector:
+  R = np.array([
+    [np.cos(alpha), -np.sin(alpha)],
+    [np.sin(alpha), np.cos(alpha)],
+  ])
+
+  return Vector(R @ v.a)
+
+
+def rotate3d(v: Vector, alpha: float, axis: int) -> Vector:
+  R = np.eye(3)
+
+  if axis == 0:
+    R[1, 1] = np.cos(alpha)
+    R[1, 2] = -np.sin(alpha)
+    R[2, 1] = np.sin(alpha)
+    R[2, 2] = np.cos(alpha)
+
+  elif axis == 1:
+    R[0, 0] = np.cos(alpha)
+    R[0, 2] = np.sin(alpha)
+    R[2, 0] = -np.sin(alpha)
+    R[2, 2] = np.cos(alpha)
+
+  elif axis == 2:
+    R[0, 0] = np.cos(alpha)
+    R[0, 1] = -np.sin(alpha)
+    R[1, 0] = np.sin(alpha)
+    R[1, 1] = np.cos(alpha)
+
+  else:
+    raise ValueError
+
+  return Vector(R @ v.a)
+
+
+def rotate_rodrigues(v: Vector, k: Vector, alpha: float) -> Vector:
+  k = k.normalize()
+
+  return v * np.cos(alpha) + (k | v) * np.sin(alpha) + k * (k @ v) * (1 - np.cos(alpha))

src/geometry/sphere.py

@@ -0,0 +1,66 @@
+from dataclasses import dataclass
+from typing import Optional
+
+import numpy as np
+
+from geometry.line_drawers import Painter
+from geometry.ray import Ray
+from geometry.rotations import rotate_rodrigues
+from geometry.surface import Surface
+from util import double_equality, minpos, solve_quadratic
+from geometry.vector import Vector
+
+
+@dataclass
+class Sphere(Surface):
+  refraction_index: float
+  r0: Vector
+  radius: float
+
+  def __and__(self, v: Vector) -> bool:
+    return double_equality(abs(self.r0 - v), self.radius)
+
+  def __contains__(self, v: Vector) -> bool:
+    return self & v or abs(self.r0 - v) < self.radius
+
+  def normal_at(self, v: Vector) -> Vector:
+    return (v - self.r0).normalize()
+
+  def intersect(self, ray: Ray) -> Optional[float]:
+    b = -2 * ((self.r0 - ray.r0) @ ray.e)
+    c = (self.r0 - ray.r0) @ (self.r0 - ray.r0) - self.radius ** 2
+
+    return minpos(solve_quadratic(1, b, c))
+
+  def draw3d(self, painter: Painter):
+    a1 = Vector.fromargs(0, 0, 1)
+    a2 = Vector.fromargs(1, 0, 0)
+    a3 = Vector.fromargs(0, 1, 0)
+
+    painter.line3d(self.r0, self.r0 + a1, lw=5, zorder=0, color='red')
+    painter.line3d(self.r0, self.r0 + a2, lw=5, zorder=0, color='green')
+    painter.line3d(self.r0, self.r0 + a3, lw=5, zorder=0, color='blue')
+
+    N1 = 8
+    N2 = 100
+
+    for lat in np.linspace(0, np.pi, N1):
+      vs = []
+      e = rotate_rodrigues(a1, a2, lat)
+
+      for long in np.linspace(0, 2 * np.pi, N2):
+        r = Ray(r0=self.r0, e=rotate_rodrigues(e, a1, long))
+        vs.append(r.prolong(self.intersect(r)))
+
+      painter.plot3d(vs, color='navy')
+
+    for long in np.linspace(0, np.pi, N1):
+      vs = []
+      e = rotate_rodrigues(a2, a1, long)
+      a3_e = rotate_rodrigues(a3, a1, long)
+
+      for lat in np.linspace(0, 2 * np.pi, N2):
+        r = Ray(r0=self.r0, e=rotate_rodrigues(e, a3_e, lat))
+        vs.append(r.prolong(self.intersect(r)))
+
+      painter.plot3d(vs, color='indigo')

src/geometry/surface.py

@@ -0,0 +1,31 @@
+from dataclasses import dataclass
+from typing import Optional
+
+from geometry.ray import Ray
+from geometry.vector import Vector
+
+
+@dataclass
+class Surface:
+  refraction_index: float
+
+  def __and__(self, v: Vector) -> bool:
+    """
+    Does this point lie on the surface itself?
+    """
+    raise NotImplementedError()
+
+  def __contains__(self, v: Vector) -> bool:
+    """
+    Is this point _inside_ the region of space, separated by this surface?
+    """
+    raise NotImplementedError()
+
+  def normal_at(self, v: Vector) -> Vector:
+    raise NotImplementedError()
+
+  def intersect(self, ray: Ray) -> Optional[float]:
+    raise NotImplementedError()
+
+  def draw2d(self, plt, **kwargs):
+    raise NotImplementedError()

src/geometry/transformed_surface.py

@@ -0,0 +1,53 @@
+from dataclasses import dataclass
+from typing import Optional
+
+import numpy as np
+
+from geometry.affine import AffineTransform
+from geometry.ray import Ray
+from geometry.surface import Surface
+from geometry.line_drawers import draw_line_3d, Painter, TransformedPainter
+from geometry.vector import Vector
+from geometry.rotations import rotate2d, rotate3d, rotate_rodrigues
+from util import double_equality, minpos, solve_quadratic
+
+
+@dataclass
+class TransformedSurface(Surface):
+  refraction_index: float
+  subsurface: Surface
+  transform: AffineTransform
+
+  def __contains__(self, v: Vector) -> bool:
+    return self.transform.apply_shift(v) in self.subsurface
+
+  def __and__(self, v: Vector) -> bool:
+    return self.subsurface & self.transform.apply_shift(v)
+
+  def normal_at(self, v: Vector) -> Vector:
+    v_ = self.transform.apply_shift(v)
+    result_ = self.subsurface.normal_at(v_)
+    result = self.transform.invert_no_shift(result_)
+    return result.normalize()
+
+  def intersect(self, ray: Ray) -> Optional[float]:
+    e_ = self.transform.apply_no_shift(ray.e)
+    length_ratio = abs(e_)
+
+    ray_ = Ray(
+      r0=self.transform.apply_shift(ray.r0),
+      e=e_.normalize(),
+    )
+
+    t_ = self.subsurface.intersect(ray_)
+    if t_ is None:
+      return None
+
+    t = t_ / length_ratio
+
+    return t
+
+  def draw3d(self, painter: Painter):
+    transformed_painter = TransformedPainter(painter, self.transform)
+
+    self.subsurface.draw3d(transformed_painter)

src/geometry/vector.py

@@ -0,0 +1,104 @@
+from __future__ import annotations
+
+import numpy as np
+
+from util import iszero
+
+
+class Vector:
+  def __init__(self, a):
+    self.a: np.ndarray = np.array(a, dtype=float)
+    self.a.setflags(write=False)
+
+    assert len(self.a.shape) == 2
+
+  @staticmethod
+  def fromiterable(iterable):
+    a = np.zeros((len(iterable), 1), dtype=float)
+
+    for i, j in enumerate(iterable):
+      a[i, 0] = j
+
+    return Vector(a)
+
+  @staticmethod
+  def fromargs(*args):
+    return Vector.fromiterable(args)
+
+  @staticmethod
+  def random():
+    return Vector.fromiterable(np.random.uniform(low=-10, high=10, size=3))
+
+  def __add__(self, other: Vector) -> Vector:
+    return Vector(self.a + other.a)
+
+  def __sub__(self, other: Vector) -> Vector:
+    return Vector(self.a - other.a)
+
+  def __matmul__(self, other: Vector) -> float:
+    """
+    Dot product
+    """
+    # noinspection PyTypeChecker
+    return np.sum(self.a * other.a)
+
+  def __or__(self, other: Vector) -> Vector:
+    """
+    Cross product
+    """
+
+    if self.a.shape != other.a.shape or self.a.shape != (3, 1):
+      raise ValueError()
+
+    return Vector.fromargs(
+      self.a[1, 0] * other.a[2, 0] - self.a[2, 0] * other.a[1, 0],
+      self.a[2, 0] * other.a[0, 0] - self.a[0, 0] * other.a[2, 0],
+      self.a[0, 0] * other.a[1, 0] - self.a[1, 0] * other.a[0, 0],
+    )
+
+  def __abs__(self) -> float:
+    """
+    Norm of self
+    """
+    return np.sqrt(self @ self)
+
+  def __eq__(self, other: Vector):
+    if self.a.shape != other.a.shape:
+      return False
+
+    return (self - other).iszero()
+
+  def __str__(self):
+    return str(self.a)
+
+  def __repr__(self):
+    return f'Vector({repr(self.a)})'
+
+  def __mul__(self, other):
+    return Vector(self.a * other)
+
+  def __rmul__(self, other):
+    return Vector(other * self.a)
+
+  def __truediv__(self, other):
+    return Vector(self.a / other)
+
+  def __mod__(self, other):
+    """
+    Is self perpendicular to other?
+    """
+    return iszero(self @ other)
+
+  def __neg__(self):
+    return -1 * self
+
+  def iszero(self):
+    return iszero(self @ self)
+
+  def normalize(self):
+    if self.iszero():
+      raise ValueError()
+    return self / abs(self)
+
+  def __hash__(self):
+    return hash(self.a)

src/rendering/_3d/render.py

@@ -0,0 +1,80 @@
+import numpy as np
+from matplotlib import pyplot as plt
+# noinspection PyUnresolvedReferences
+from mpl_toolkits.mplot3d import Axes3D
+
+from geometry.ellipse import Ellipse
+from geometry.line_drawers import draw_line_3d
+from geometry.ray import Ray
+from geometry.rotations import rotate_rodrigues
+from geometry.vector import Vector
+from rendering import raytracing
+from rendering.raytracing import RaytracingSettings
+
+figure = plt.figure(figsize=(8, 8), dpi=200)
+axis = figure.add_subplot(111, projection='3d')
+axis.view_init(elev=10, azim=-80)
+
+surfaces = []
+
+ellipse = Ellipse(
+  refraction_index=1.5,
+  f1=Vector.fromargs(-2, 2, 1),
+  f2=Vector.fromargs(1, 0, -2),
+  distance=5.8,
+)
+surfaces.append(ellipse)
+ellipse.draw3d(axis)
+
+axis.set_xlim(-4, 4)
+axis.set_ylim(-4, 4)
+axis.set_zlim(-4, 4)
+
+axis.set_xlabel('x')
+axis.set_ylabel('y')
+axis.set_zlabel('z')
+
+axis.grid(True)
+
+settings = RaytracingSettings(
+  STARTING_LW=2,
+  STARTING_HUE=0,
+
+  REFLECT_HUE_SHIFT=0,
+  REFLECT_LW_DECAY=0.7,
+
+  REFRACT_HUE_SHIFT=1 / 8,
+  REFRACT_LW_DECAY=0.8,
+
+  MAX_DEPTH=3
+)
+
+for alpha in np.linspace(-np.pi / 32, np.pi / 32, 4):
+  for beta in np.linspace(-np.pi / 32, np.pi / 32, 4):
+    e = rotate_rodrigues(
+      Vector.fromargs(-1, 0, 0),
+      Vector.fromargs(0, 1, 0),
+      alpha,
+    )
+    e = rotate_rodrigues(
+      e,
+      Vector.fromargs(0, 0, 1),
+      beta,
+    )
+
+    raytracing.raytrace(
+      Ray(
+        r0=Vector.fromargs(4, 0, 0),
+        e=e
+      ),
+      settings,
+      axis,
+      surfaces,
+      draw_line_3d
+    )
+
+figure.show()
+
+for i, azim in enumerate(np.linspace(0, 360, 16)):
+  axis.view_init(elev=10, azim=azim)
+  figure.savefig(f'movie/frame{i:03d}.png')

src/rendering/_3d/render_transformed.py

@@ -0,0 +1,111 @@
+import numpy as np
+from matplotlib import pyplot as plt
+# noinspection PyUnresolvedReferences
+from mpl_toolkits.mplot3d import Axes3D
+
+from geometry.affine import AffineTransform
+from geometry.line_drawers import draw_line_3d, Painter
+from geometry.ray import Ray
+from geometry.rotations import rotate_rodrigues
+from geometry.sphere import Sphere
+from geometry.transformed_surface import TransformedSurface
+from geometry.vector import Vector
+from rendering import raytracing
+from rendering.raytracing import RaytracingSettings
+
+figure = plt.figure(figsize=(8, 8), dpi=200)
+axis = figure.add_subplot(111, projection='3d')
+axis.view_init(elev=10, azim=-80)
+
+painter = Painter(axis)
+
+surfaces = []
+
+# sphere = Sphere(
+#   refraction_index=1.5,
+#   r0=Vector.fromargs(0, 0, 1),
+#   radius=2,
+# )
+# surfaces.append(sphere)
+# sphere.draw3d(painter)
+
+np.random.seed(5)
+transform_matrix = np.random.normal(loc=1, scale=2, size=(3, 3))
+transform_shift = np.random.normal(loc=0, scale=2, size=(3, 1))
+
+transform = AffineTransform(
+  # A=np.array([
+  #   [2, 0, 0],
+  #   [0, 0.5, 0],
+  #   [0, 0, 1],
+  # ]),
+  A=transform_matrix,
+  # b=Vector.fromargs(0, 0, 0),
+  b=Vector(transform_shift),
+)
+sphere = Sphere(
+  refraction_index=1.5,
+  r0=Vector.fromargs(0, 0, 1),
+  radius=1,
+)
+transformed_sphere = TransformedSurface(
+  refraction_index=sphere.refraction_index,
+  subsurface=sphere,
+  transform=transform
+)
+# sphere.draw3d(painter)
+transformed_sphere.draw3d(painter)
+surfaces.append(transformed_sphere)
+
+axis.set_xlim(-4, 4)
+axis.set_ylim(-4, 4)
+axis.set_zlim(-4, 4)
+
+axis.set_xlabel('x')
+axis.set_ylabel('y')
+axis.set_zlabel('z')
+
+axis.grid(True)
+
+settings = RaytracingSettings(
+  STARTING_LW=2,
+  STARTING_HUE=0,
+
+  REFLECT_HUE_SHIFT=0,
+  REFLECT_LW_DECAY=0.7,
+
+  REFRACT_HUE_SHIFT=1 / 8,
+  REFRACT_LW_DECAY=0.8,
+
+  MAX_DEPTH=3
+)
+
+for alpha in np.linspace(-np.pi / 32, np.pi / 32, 4):
+  for beta in np.linspace(-np.pi / 32, np.pi / 32, 4):
+    e = rotate_rodrigues(
+      Vector.fromargs(-1, 0, 0),
+      Vector.fromargs(0, 1, 0),
+      alpha,
+    )
+    e = rotate_rodrigues(
+      e,
+      Vector.fromargs(0, 0, 1),
+      beta,
+    )
+
+    raytracing.raytrace(
+      Ray(
+        r0=Vector.fromargs(4, 0, 0),
+        e=e
+      ),
+      settings,
+      axis,
+      surfaces,
+      draw_line_3d
+    )
+
+figure.show()
+
+for i, azim in enumerate(np.linspace(0, 360, 16)):
+  axis.view_init(elev=10, azim=azim)
+  figure.savefig(f'movie/frame{i:03d}.png')

src/rendering/raytracing.py

@@ -0,0 +1,110 @@
+import dataclasses
+from typing import List
+
+from matplotlib.colors import hsv_to_rgb
+
+from geometry.ray import Ray
+from geometry.reflections import reflect, refract
+from geometry.surface import Surface
+from util import iszero
+
+
+@dataclasses.dataclass
+class RaytracingSettings:
+  STARTING_LW: float
+  REFLECT_LW_DECAY: float
+  REFRACT_LW_DECAY: float
+
+  STARTING_HUE: float
+  REFLECT_HUE_SHIFT: float
+  REFRACT_HUE_SHIFT: float
+
+  MAX_DEPTH: int
+
+  lw: float = None
+  hue: float = None
+  depth: float = None
+  n: float = 1
+
+  def __post_init__(self):
+    if self.lw is None:
+      self.lw = self.STARTING_LW
+
+    if self.hue is None:
+      self.hue = self.STARTING_HUE
+
+    if self.depth is None:
+      self.depth = self.MAX_DEPTH
+
+  def copy(self):
+    return dataclasses.replace(self)
+
+  def reflect(self):
+    return dataclasses.replace(
+      self,
+      lw=self.lw * self.REFLECT_LW_DECAY,
+      hue=(self.hue + self.REFLECT_HUE_SHIFT) % 1,
+      depth=self.depth - 1,
+    )
+
+  def refract(self, new_n):
+    return dataclasses.replace(
+      self,
+      lw=self.lw * self.REFRACT_LW_DECAY,
+      hue=(self.hue + self.REFRACT_HUE_SHIFT) % 1,
+      depth=self.depth - 1,
+      n=new_n,
+    )
+
+
+def raytrace(ray: Ray,
+             settings: RaytracingSettings,
+             axis,
+             surfaces: List[Surface],
+             line_drawer):
+  if settings.depth == 0:
+    return
+
+  min_t = None
+  surface = None
+
+  for cur_surface in surfaces:
+    t = cur_surface.intersect(ray)
+    if t is None or iszero(t):
+      continue
+
+    if min_t is None or t < min_t:
+      min_t = t
+      surface = cur_surface
+
+  if min_t is None:
+    line_drawer(axis, ray.r0, ray.prolong(1),
+                color=hsv_to_rgb([settings.hue, settings.lw / settings.STARTING_LW, 1]),
+                zorder=settings.depth + 5,
+                linestyle='--',
+                lw=settings.lw)
+    return
+
+  intersection = ray.prolong(min_t)
+  line_drawer(axis, ray.r0, intersection,
+              color=hsv_to_rgb([settings.hue, settings.lw / settings.STARTING_LW, 1]),
+              zorder=settings.depth + 5,
+              linestyle='-',
+              lw=settings.lw)
+
+  if ray.r0 not in surface:
+    # moves into
+    normal = surface.normal_at(intersection)
+    new_n = settings.n * surface.refraction_index
+  else:
+    # moves outside
+    normal = -surface.normal_at(intersection)
+    new_n = settings.n / surface.refraction_index
+
+  # reflection
+  reflected = reflect(ray.e, normal)
+  raytrace(Ray(intersection, reflected), settings.reflect(), axis, surfaces, line_drawer)
+
+  # refraction
+  refracted = refract(ray.e, normal, settings.n, new_n)
+  raytrace(Ray(intersection, refracted), settings.refract(new_n), axis, surfaces, line_drawer)

src/rendering/render.py

@@ -0,0 +1,111 @@
+from typing import Set
+
+import numpy as np
+from matplotlib import pyplot as plt
+
+from plane import Plane
+from ray import Ray
+from reflections import reflect, refract
+from sphere import Sphere
+from vector import Vector
+from util import iszero
+
+
+def plot_v(v1: Vector, v2: Vector, *args, **kwargs):
+  plt.plot(
+    [v1.a[0, 0], v2.a[0, 0]],
+    [v1.a[1, 0], v2.a[1, 0]],
+    *args, **kwargs)
+
+
+intersectables = []
+f = plt.figure(figsize=(6, 6))
+ax = plt.gca()
+
+# declare a plane
+plane = Plane(
+  r0=Vector.fromargs(1, 0),
+  n=Vector.fromargs(-2, 1),
+)
+intersectables.append((plane, 1.5))
+plt.plot([-1, 3], [-4, 4], label='плоскость')
+
+# declare a sphere
+sphere = Sphere(
+  r0=Vector.fromargs(0, 2),
+  radius=0.8,
+)
+intersectables.append((sphere, 2))
+circle = plt.Circle(
+  (sphere.r0.a[0, 0], sphere.r0.a[1, 0]),
+  sphere.radius, color='blue', fill=False)
+ax.add_artist(circle)
+
+
+def raytrace(
+    ray: Ray,
+    current_n: float,
+    already_intersected: Set,
+    depth: int):
+  if depth == 0:
+    return
+
+  min_t = None
+  surface = None
+  surface_n = None
+
+  for intersectable, intersectable_n in intersectables:
+    t = intersectable.intersect(ray)
+    if t is None or iszero(t):
+      continue
+
+    if min_t is None or t < min_t:
+      min_t = t
+      surface = intersectable
+      surface_n = intersectable_n
+
+  if min_t is None:
+    plot_v(ray.r0, ray.prolong(5), color='red')
+    return
+
+  intersection = ray.prolong(min_t)
+  plot_v(ray.r0, intersection, color='orange')
+
+  normal = surface.normal_at(intersection)
+  new_n = surface_n
+  new_already_intersected = set(already_intersected)
+
+  if id(surface) not in already_intersected:
+    # moves_into
+    new_already_intersected.add(id(surface))
+  else:
+    # moves outside
+    new_already_intersected.remove(id(surface))
+    normal = -1 * normal
+    new_n = 1
+
+  # reflection
+  reflected = reflect(ray.e, normal)
+  raytrace(Ray(intersection, reflected), current_n, already_intersected, depth - 1)
+
+  # refraction
+  refracted = refract(ray.e, normal, current_n, new_n)
+  raytrace(Ray(intersection, refracted), new_n, new_already_intersected, depth - 1)
+
+
+raytrace(Ray(
+  r0=Vector.fromargs(-4, 2),
+  e=Vector.fromargs(2, -0.5)
+), 1, set(), 10)
+
+# for i in range(30):
+#   trace_ray(Ray(
+#     r0=Vector.fromargs(-4, -2),
+#     e=Vector.fromargs(np.cos(i / 20 / np.pi), np.sin(i / 20 / np.pi)),
+#   ))
+
+plt.xlim(-5, 5)
+plt.ylim(-5, 5)
+# plt.legend()
+plt.grid(True)
+plt.show()

src/rendering/render2.py

@@ -0,0 +1,125 @@
+from typing import Set
+
+import numpy as np
+from matplotlib import pyplot as plt
+from matplotlib.colors import hsv_to_rgb
+
+from plane import Plane
+from ray import Ray
+from reflections import reflect, refract
+from sphere import Sphere
+from vector import Vector
+from util import iszero
+
+
+def draw_line(v1: Vector, v2: Vector, **kwargs):
+  plt.plot(
+    [v1.a[0, 0], v2.a[0, 0]],
+    [v1.a[1, 0], v2.a[1, 0]],
+    **kwargs)
+
+
+def draw_point(v: Vector, **kwargs):
+  plt.gca().add_artist(
+    plt.Circle(
+      (v.a[0, 0], v.a[1, 0]), 0.1,
+      fill=False, **kwargs))
+
+
+f = plt.figure(figsize=(8, 8))
+
+surfaces = []
+
+np.random.seed(5)
+
+for i in range(5):
+  sphere = Sphere(
+    refraction_index=3,
+    r0=Vector.fromiterable(np.random.uniform(-2, 2, size=2)),
+    radius=np.random.uniform(0.5, 2),
+  )
+  surfaces.append(sphere)
+
+  plt.gca().add_artist(plt.Circle(
+    (sphere.r0.a[0, 0], sphere.r0.a[1, 0]),
+    sphere.radius, color='blue', fill=False, zorder=200))
+  plt.gca().add_artist(plt.Circle(
+    (sphere.r0.a[0, 0], sphere.r0.a[1, 0]),
+    sphere.radius, color='blue', fill=True, zorder=0, alpha=0.05))
+
+
+def raytrace(
+    ray: Ray,
+    hue: float,
+    lw: float,
+    current_n: float,
+    depth: int):
+  if depth == 0:
+    return
+
+  min_t = None
+  surface = None
+
+  for cur_surface in surfaces:
+    t = cur_surface.intersect(ray)
+    if t is None or iszero(t):
+      continue
+
+    if min_t is None or t < min_t:
+      min_t = t
+      surface = cur_surface
+
+  if min_t is None:
+    draw_line(ray.r0, ray.prolong(5), color=hsv_to_rgb([hue, 1, 1]), zorder=depth + 5, lw=lw)
+    return
+
+  intersection = ray.prolong(min_t)
+  draw_line(ray.r0, intersection, color=hsv_to_rgb([hue, 1, 1]), zorder=depth + 5, lw=lw)
+
+  # normal = surface.normal_at(intersection)
+  # new_n = current_n *
+
+  if ray.r0 not in surface:
+    # moves into
+    normal = surface.normal_at(intersection)
+    new_n = current_n * surface.refraction_index
+  else:
+    # moves outside
+    normal = -1 * surface.normal_at(intersection)
+    new_n = current_n / surface.refraction_index
+
+  # reflection
+  reflected = reflect(ray.e, normal)
+  raytrace(
+    Ray(intersection, reflected),
+    hue=hue,
+    lw=lw * 1,
+    current_n=current_n,
+    depth=depth - 1)
+
+  # refraction
+  refracted = refract(ray.e, normal, current_n, new_n)
+  raytrace(
+    Ray(intersection, refracted),
+    hue=(hue + 0.15) % 1,
+    lw=lw * 0.7,
+    current_n=new_n,
+    depth=depth - 1)
+
+
+raytrace(Ray(
+  r0=Vector.fromargs(-3.5, 0),
+  e=Vector.fromargs(2, -0.5)
+), hue=0, lw=5, current_n=1, depth=8)
+
+# for i in range(30):
+#   trace_ray(Ray(
+#     r0=Vector.fromargs(-4, -2),
+#     e=Vector.fromargs(np.cos(i / 20 / np.pi), np.sin(i / 20 / np.pi)),
+#   ))
+
+plt.xlim(-4, 4)
+plt.ylim(-4, 4)
+# plt.legend()
+plt.grid(True, zorder=1)
+plt.show()

src/rendering/render3.py

@@ -0,0 +1,131 @@
+import numpy as np
+from matplotlib import pyplot as plt
+from matplotlib.colors import hsv_to_rgb
+
+from ellipse import Ellipse
+from ray import Ray
+from reflections import reflect, refract
+from util import iszero
+from vector import Vector
+
+MAX_DEPTH = 10
+
+STARTING_LW = 5
+
+REFRACT_HUE_SHIFT = 1 / 8
+
+LW_DECAY = 0.8
+
+
+def draw_line(v1: Vector, v2: Vector, **kwargs):
+  plt.plot(
+    [v1.a[0, 0], v2.a[0, 0]],
+    [v1.a[1, 0], v2.a[1, 0]],
+    **kwargs)
+
+
+def draw_point(v: Vector, **kwargs):
+  plt.gca().add_artist(
+    plt.Circle(
+      (v.a[0, 0], v.a[1, 0]), 0.1,
+      fill=False, **kwargs))
+
+
+f = plt.figure(figsize=(8, 8))
+
+surfaces = []
+
+np.random.seed(2)
+
+# ellipse = Ellipse(
+#   refraction_index=1.3,
+#   f1=Vector.fromargs(1, 0),
+#   f2=Vector.fromargs(-1, 0),
+#   distance=3,
+# )
+# ellipse.draw2d(plt)
+# surfaces.append(ellipse)
+
+for i in range(5):
+  f1 = Vector.fromiterable(np.random.uniform(-3, 3, size=2))
+  f2 = Vector.fromiterable(np.random.uniform(-3, 3, size=2))
+  distance = abs(f1 - f2) + np.random.uniform(0, 2)
+  ellipse = Ellipse(
+    refraction_index=2,
+    f1=f1,
+    f2=f2,
+    distance=distance,
+  )
+  ellipse.draw2d(plt, color='blue')
+  surfaces.append(ellipse)
+
+
+def raytrace(
+    ray: Ray,
+    hue: float,
+    lw: float,
+    current_n: float,
+    depth: int):
+  if depth == 0:
+    return
+
+  min_t = None
+  surface = None
+
+  for cur_surface in surfaces:
+    t = cur_surface.intersect(ray)
+    if t is None or iszero(t):
+      continue
+
+    if min_t is None or t < min_t:
+      min_t = t
+      surface = cur_surface
+
+  if min_t is None:
+    draw_line(ray.r0, ray.prolong(5), color=hsv_to_rgb([hue, 1, 1]), zorder=depth + 5, lw=lw)
+    return
+
+  intersection = ray.prolong(min_t)
+  draw_line(ray.r0, intersection, color=hsv_to_rgb([hue, 1, 1]), zorder=depth + 5, lw=lw)
+
+  # normal = surface.normal_at(intersection)
+  # new_n = current_n *
+
+  if ray.r0 not in surface:
+    # moves into
+    normal = surface.normal_at(intersection)
+    new_n = current_n * surface.refraction_index
+  else:
+    # moves outside
+    normal = -1 * surface.normal_at(intersection)
+    new_n = current_n / surface.refraction_index
+
+  # reflection
+  reflected = reflect(ray.e, normal)
+  raytrace(
+    Ray(intersection, reflected),
+    hue=hue,
+    lw=lw * LW_DECAY,
+    current_n=current_n,
+    depth=depth - 1)
+
+  # refraction
+  refracted = refract(ray.e, normal, current_n, new_n)
+  raytrace(
+    Ray(intersection, refracted),
+    hue=(hue + REFRACT_HUE_SHIFT) % 1,
+    lw=lw * LW_DECAY,
+    current_n=new_n,
+    depth=depth - 1)
+
+
+raytrace(Ray(
+  r0=Vector.fromargs(-3, 3),
+  e=Vector.fromargs(1, -1)
+), hue=0, lw=STARTING_LW, current_n=1, depth=MAX_DEPTH)
+
+plt.xlim(-4, 4)
+plt.ylim(-4, 4)
+# plt.legend()
+plt.grid(True, zorder=1)
+plt.show()

src/rendering/render4.py

@@ -0,0 +1,136 @@
+import numpy as np
+from matplotlib import pyplot as plt
+from matplotlib.colors import hsv_to_rgb
+
+from geometry.ellipse import Ellipse
+from geometry.ray import Ray
+from geometry.reflections import reflect, refract
+from geometry.vector import Vector
+from geometry.rotations import rotate2d
+from util import iszero
+
+MAX_DEPTH = 2
+
+STARTING_LW = 3
+
+REFRACT_HUE_SHIFT = 1 / 8
+
+LW_DECAY = 0.7
+
+
+def draw_line(v1: Vector, v2: Vector, **kwargs):
+  plt.plot(
+    [v1.a[0, 0], v2.a[0, 0]],
+    [v1.a[1, 0], v2.a[1, 0]],
+    **kwargs)
+
+
+def draw_point(v: Vector, **kwargs):
+  plt.gca().add_artist(
+    plt.Circle(
+      (v.a[0, 0], v.a[1, 0]), 0.1,
+      fill=False, **kwargs))
+
+
+f = plt.figure(figsize=(8, 8), dpi=200)
+
+surfaces = []
+
+np.random.seed(2)
+
+
+# sphere = Sphere(
+#   refraction_index=1.3,
+#   r0=Vector.fromargs(-1, 2),
+#   radius=1,
+# )
+# # sphere.draw2d(plt)
+# surfaces.append(sphere)
+
+
+def raytrace(
+    ray: Ray,
+    hue: float,
+    lw: float,
+    current_n: float,
+    depth: int):
+  if depth == 0:
+    return
+
+  min_t = None
+  surface = None
+
+  for cur_surface in surfaces:
+    t = cur_surface.intersect(ray)
+    if t is None or iszero(t):
+      continue
+
+    if min_t is None or t < min_t:
+      min_t = t
+      surface = cur_surface
+
+  if min_t is None:
+    draw_line(ray.r0, ray.prolong(10),
+              color=hsv_to_rgb([hue, lw / STARTING_LW, 1]),
+              zorder=depth + 5, lw=lw)
+    return
+
+  intersection = ray.prolong(min_t)
+  draw_line(ray.r0, intersection,
+            color=hsv_to_rgb([hue, lw / STARTING_LW, 1]),
+            zorder=depth + 5, lw=lw)
+
+  # normal = surface.normal_at(intersection)
+  # new_n = current_n *
+
+  if ray.r0 not in surface:
+    # moves into
+    normal = surface.normal_at(intersection)
+    new_n = current_n * surface.refraction_index
+  else:
+    # moves outside
+    normal = -1 * surface.normal_at(intersection)
+    new_n = current_n / surface.refraction_index
+
+  # reflection
+  reflected = reflect(ray.e, normal)
+  raytrace(
+    Ray(intersection, reflected),
+    hue=hue,
+    lw=lw * LW_DECAY,
+    current_n=current_n,
+    depth=depth - 1)
+
+  # refraction
+  refracted = refract(ray.e, normal, current_n, new_n)
+  raytrace(
+    Ray(intersection, refracted),
+    hue=(hue + REFRACT_HUE_SHIFT) % 1,
+    lw=lw * LW_DECAY,
+    current_n=new_n,
+    depth=depth - 1)
+
+
+ellipse = Ellipse(
+  refraction_index=1,
+  f1=Vector.fromargs(0, -2),
+  f2=Vector.fromargs(2, 2),
+  distance=5,
+)
+ellipse.draw2d(plt, color='blue')
+surfaces.append(ellipse)
+
+r0 = Vector.fromargs(-3, 0)
+e = Vector.fromargs(1, 0)
+
+for alpha in np.linspace(-np.pi / 4, np.pi / 4, 31):
+  raytrace(Ray(
+    r0=r0,
+    e=rotate2d(e, alpha)
+  ), hue=0, lw=STARTING_LW, current_n=1, depth=MAX_DEPTH)
+
+plt.xlim(-4, 4)
+plt.ylim(-4, 4)
+# plt.legend()
+plt.grid(True, zorder=1)
+plt.show()

src/rendering/render5.py

@@ -0,0 +1,142 @@
+import numpy as np
+from matplotlib import pyplot as plt
+from matplotlib.colors import hsv_to_rgb
+
+from ellipse import Ellipse
+from ray import Ray
+from reflections import reflect, refract
+from util import iszero
+from vector import Vector, rotate2d
+
+MAX_DEPTH = 6
+
+STARTING_LW = 5
+
+REFRACT_HUE_SHIFT = 1 / 8
+
+LW_DECAY = 0.7
+
+
+def draw_line(v1: Vector, v2: Vector, **kwargs):
+  plt.plot(
+    [v1.a[0, 0], v2.a[0, 0]],
+    [v1.a[1, 0], v2.a[1, 0]],
+    **kwargs)
+
+
+def draw_point(v: Vector, **kwargs):
+  plt.gca().add_artist(
+    plt.Circle(
+      (v.a[0, 0], v.a[1, 0]), 0.1,
+      fill=False, **kwargs))
+
+
+f = plt.figure(figsize=(8, 8), dpi=200)
+
+surfaces = []
+
+np.random.seed(2)
+
+
+# ellipse = Ellipse(
+#   refraction_index=1.3,
+#   f1=Vector.fromargs(1, 0),
+#   f2=Vector.fromargs(-1, 0),
+#   distance=3,
+# )
+# ellipse.draw2d(plt)
+# surfaces.append(ellipse)
+
+
+def raytrace(
+    ray: Ray,
+    hue: float,
+    lw: float,
+    current_n: float,
+    depth: int):
+  if depth == 0:
+    return
+
+  min_t = None
+  surface = None
+
+  for cur_surface in surfaces:
+    t = cur_surface.intersect(ray)
+    if t is None or iszero(t):
+      continue
+
+    if min_t is None or t < min_t:
+      min_t = t
+      surface = cur_surface
+
+  if min_t is None:
+    draw_line(ray.r0, ray.prolong(10),
+              color=hsv_to_rgb([hue, lw / STARTING_LW, 1]),
+              zorder=depth + 5, lw=lw)
+    return
+
+  intersection = ray.prolong(min_t)
+  draw_line(ray.r0, intersection,
+            color=hsv_to_rgb([hue, lw / STARTING_LW, 1]),
+            zorder=depth + 5, lw=lw)
+
+  if ray.r0 not in surface:
+    # moves into
+    normal = surface.normal_at(intersection)
+    new_n = current_n * surface.refraction_index
+  else:
+    # moves outside
+    normal = -1 * surface.normal_at(intersection)
+    new_n = current_n / surface.refraction_index
+
+  # reflection
+  reflected = reflect(ray.e, normal)
+  raytrace(
+    Ray(intersection, reflected),
+    hue=hue,
+    lw=lw * LW_DECAY,
+    current_n=current_n,
+    depth=depth - 1)
+
+  # refraction
+  refracted = refract(ray.e, normal, current_n, new_n)
+  raytrace(
+    Ray(intersection, refracted),
+    hue=(hue + REFRACT_HUE_SHIFT) % 1,
+    lw=lw * LW_DECAY,
+    current_n=new_n,
+    depth=depth - 1)
+
+
+ellipse1 = Ellipse(
+  refraction_index=3,
+  f1=Vector.fromargs(2, -2),
+  f2=Vector.fromargs(3, 2),
+  distance=5,
+)
+ellipse1.draw2d(plt, color='blue')
+surfaces.append(ellipse1)
+
+ellipse2 = Ellipse(
+  refraction_index=3,
+  f1=Vector.fromargs(-1, -2),
+  f2=Vector.fromargs(-1, 2),
+  distance=4.2,
+)
+ellipse2.draw2d(plt, color='blue')
+surfaces.append(ellipse2)
+
+r0 = Vector.fromargs(-3, 0)
+e = Vector.fromargs(1, 0)
+
+for alpha in np.linspace(-np.pi / 4, np.pi / 4, 9):
+  raytrace(Ray(
+    r0=r0,
+    e=rotate2d(e, alpha)
+  ), hue=0, lw=STARTING_LW, current_n=1, depth=MAX_DEPTH)
+
+plt.xlim(-4, 4)
+plt.ylim(-4, 4)
+# plt.legend()
+plt.grid(True, zorder=1)
+plt.show()

src/test.py

@@ -0,0 +1,18 @@
+import numpy as np
+
+from plane import Plane
+from vector import Vector
+
+if __name__ == '__main__':
+  a = Vector.fromargs(1, 2, 3)
+  b = Vector.fromiterable(np.eye(3)[:, 0])
+  print(a @ b)
+  print(abs(a))
+  print(a / 2)
+
+  p = Plane(
+    Vector.fromargs(0, 0, 0),
+    Vector.fromargs(0, 1, 0),
+  )
+
+  print(p & Vector.fromargs(1, 1, 0))

src/util.py

@@ -0,0 +1,46 @@
+from __future__ import annotations
+
+from typing import List, Optional
+
+import numpy as np
+
+EPS = 1e-6
+
+
+def double_equality(a, b) -> bool:
+  return abs(a - b) < EPS
+
+
+def iszero(x) -> bool:
+  return double_equality(np.sum(x), 0)
+
+
+def minpos_criterion(x: float) -> bool:
+  return x is not None and not iszero(x) and x > 0
+
+
+def minpos(args) -> Optional[float]:
+  args = list(filter(minpos_criterion, args))
+
+  if not args:
+    return None
+
+  return min(args)
+
+
+def solve_quadratic(a, b, c) -> List[float]:
+  if iszero(a):
+    if iszero(b):
+      return []
+
+    return [-c / b]
+
+  D = b ** 2 - 4 * a * c
+
+  if iszero(D):
+    return [-b / 2 / a]
+
+  if D < 0:
+    return []
+
+  return [(-b - np.sqrt(D)) / 2 / a, (-b + np.sqrt(D)) / 2 / a]

tests/test_intersections.py

@@ -0,0 +1,30 @@
+import numpy as np
+
+from intersections import intersect_plane_ray
+from plane import Plane
+from ray import Ray
+from vector import Vector
+
+
+def test_intersections():
+  r = Ray(
+    Vector.fromargs(0, 0, 1),
+    Vector.fromargs(-1, -1, -1),
+  )
+  p = Plane(
+    Vector.fromargs(0, 0, 0),
+    Vector.fromargs(0, 0, 1),
+  )
+
+  assert intersect_plane_ray(p, r) == Vector.fromargs(-1, -1, 0)
+
+  np.random.seed(0)
+  for i in range(10):
+    p = Plane(Vector.random(), Vector.random())
+    r = Ray(Vector.random(), Vector.random())
+
+    intersection = intersect_plane_ray(p, r)
+
+    print(p, r, intersection)
+
+    assert p & intersection

tests/test_plane.py

@@ -0,0 +1,37 @@
+import numpy as np
+
+from plane import Plane
+from ray import Ray
+from vector import Vector
+
+
+def test_plane():
+  p = Plane(
+    Vector.fromargs(0, 0, 0),
+    Vector.fromargs(0, 0, 1),
+  )
+  assert (p & Vector.fromargs(0, 0, 0))
+  assert (p & Vector.fromargs(1, 1, 0))
+  assert not (p & Vector.fromargs(0, 0, -1))
+
+  p = Plane(
+    Vector.fromargs(1, 1),
+    Vector.fromargs(1, 1),
+  )
+
+  assert not (p & Vector.fromargs(0, 0))
+  assert (p & Vector.fromargs(2, 0))
+
+
+def test_reflection():
+  p = Plane(
+    Vector.fromargs(0, 0),
+    Vector.fromargs(0, 1),
+  )
+
+  r = Ray(
+    Vector.fromargs(-1, 1),
+    Vector.fromargs(1, -1),
+  )
+
+  assert p.intersect(r) == np.sqrt(2)

tests/test_sphere.py

@@ -0,0 +1,22 @@
+import numpy as np
+
+from ray import Ray
+from sphere import Sphere
+from vector import Vector
+
+
+def test_intersect():
+  sphere = Sphere(
+    r0=Vector.fromargs(0, 0),
+    radius=1,
+  )
+
+  assert sphere.intersect(Ray(
+    r0=Vector.fromargs(-2, 0),
+    e=Vector.fromargs(1, 0),
+  )) == 1
+
+  assert sphere.intersect(Ray(
+    r0=Vector.fromargs(-1, -1),
+    e=Vector.fromargs(1, 1),
+  )) == np.sqrt(2) - 1

tests/test_vector.py

@@ -0,0 +1,7 @@
+from vector import Vector
+
+
+def test_vector():
+  v = Vector.fromargs(1, 2)
+  print(v)
+  print(-v)